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CALCULATION OF THE FLOW FIELD IN SUPERSONIC MIXED-COMPRESSION INLETS 
AT ANGLE OF ATTACK USING THE THREE-DIMENSIONAL METHOD OF 


CHARACTERISTICS WITH DISCRETE SHOCK WAVE FITTING* 


Joseph Vadyak and Joe D. Hoffman 
School of Mechanical Engineering 
Purdue University, West Lafayette, Indiana 47907 


SUMMARY 


An analysis has been developed for calculating the flow field in supersonic 
mixed-compression aircraft inlets at angle of attack using the three-dimension- 
al method of characteristics with discrete shock wave fitting. This report 
describes the details of the analysis and presents some computational results. 

The gas dynamic model is based on the assumptions of steady continuum flow, 
negligible body forces, a simple system in thermodynamic equilibrium, no mass 
diffusion, negligible -adiative heat transfer, no internal heat generation 
other than viscous dissipation, and viscous and thermal diffusion effects of 
secondary importance. The viscous and thermal diffusion terms are treated as 
forcing functions, or correction terms, in the method of characteristics 
scheme. Pressure and density are specified as the primary thermodynamic pro- 
perties, and the temperature, speed of sound, viscosity, and thermal conduc- 
tivity are expressed in terms of pressure and density. 

The s'^tem of governing equations is hyperbolic when the flow is supersonic. 
The equations for the characteristic surfaces and the compatibility equations 
applicable along these surfaces are derived. The characteristic surfaces are 
the stream surfaces, which are surfaces composed of streamlines, and the wave 
surfaces, which are surfaces tangent to a Mach conoid. The compatibility 
equations are expressed as directional derivatives along streamlines and bi- 
characteristics, which are the lines of tangency between a wave surface and a 
Mach conoid. The numerical integration procedure devised by D.S. Butler was 
employed to develop a numerical integration algorithm that is second-order 
accurate, explicit, and does not violate the domain of dependence of the dif- 
ferential equations. 

The bow shock wave surrounding the forebody portion of the centerbody and 
the internal shock wave system inside of the inlet are determined by discrete 
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shock wave fitting. The continuous flow field between shock waves is deter- 
mined by the method of characteristics numerical integration procedure, and the 
flow properties across the shock waves are determined by the application of the 
Hugoniot jump conditions. 

Unit processes were developed for interior field points, solid boundary 
points, field-shock wave points, and solid boundary -shock wave points. An 
inverse marching scheme is employed in which the solution is obtained on planes 
perpendicular to the axis of the centerbody and the cowl. The distance between 
successive solution planes is determined by the Courant-Friedrlchs-Lewy stabil- 
ity criterion. Although the numerical integration procedure developed herein 
is capable of analyzing three-dimensional flows in three-dimensional geometries, 
only ax i symmetric geometries at angle of attack were considered in the present 
investigation. 

Selected computational results are presented for three categories of flow 
fields: external flow about the forebody, continuous internal flow, and in- 

ternal flow in which the discrete internal shock wave system is computed. Both 
axisymmetric flow results and three-dimensional flow results are presented. For 
the internal flow field in which the shock waves have been fitted, some com- 
parisons with experimental data are presented. Results of the present analysis 
are compared with those obtained by the two-dimensional method of characteristics 
for axisyirmetric flows, and by a three-dimensional fixed grid finite difference 
shock capturing method. 


The computational results support the following conclusions. The external 
flow field about a forebody can be accurately calculated if a bow shock wave 
of reasonable strength exists. For axisymmetric flows, the solution agrees 
well with results obtained by the two-dimensional method of characteristics. 
Except in regions of strong viscous interaction and boundary layer removal, 
the results of the present analysis agree well with experimental data. Good 
agreement is obtained between the present analysis and a finite difference 
shock capturing method. The present analysis, however, which discretely fits 
shock waves, provides better resolution of the shock waves. 
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SECTION I 


INTRODUCTION 


1 . GENERAL 

The purpose of this investigation was to develop a method for 
calculating the flow field in a supersonic mixed-compression aircraft 
inlet at nonzero angle of attack. A typical supersonic mixed -compres- 
sion aircraft inlet is illustrated in Figure 1. Compression takes 
place both in the external flow about the forebody and in the internal 
flow inside the annulus. The free-stream velocity is supersonic, 
hence, a bow shock wave is generated at the forebody tip as shown. 

The internal shock wave emanates from the cowl lip. That shock wave 
makes a number of reflections with the centerbody and cowl before 
terminating in the divergence downstream of the geometric throat of 
the annulus. The flow is subsonic downstream of that location. 

A major objective in the design of any aircraft inlet is to achieve 
maximum flow compression with a minimum reduction in stagnation pres- 
sure. Moreover, since an adverse pressure gradient exists, suitable 
control of the boundary layer is a major design consideration. This 
is especially true for an inlet such as that illustrated in Figure 1, 
since a number of oblique shock wave- boundary layer interactions occur. 
In a mixed-compression inlet, it is not unusual to remove 10 percent 
or more of the cowl lip mass flow rate by boundary layer bleed to 
control separation of the boundary layer. 
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FIGURE I. MIXED-COMPRESSION AIRCRAFT INLET 
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The inlet illustrated in Figure 1 is axi symmetric. At zero 
incidence, the flow field is axi symmetric and can be computed using a 
two-dimensional method. However, at nonzero angle of attack, cross 
flow develops, and computation of the flow field requires using a three- 
dimensional algorithm. 

2. METHODS OF SOLUTION 

The equations of motion for steady three-dimensional supersonic 
flow may be classified as a system of hyperbolic quasi-linear partial 
differential equations of first order. Exact solutions can be found 
in only a few cases. Consequently, most solutions are obtained by 
employing numerical techniques. The two most widely used numerical 
methods are: 

1. method of finite differences 

2. method of characteristics 

The method of finite differences replaces the derivatives in the 
system of original differential equations with simple differences. 

The system of difference equations is then solved to obtain the solu- 
tion. Finite difference methods may be further classified into those 
methods which do and those methods which do not contain artificial 
viscosity terms. The artificial viscosity terms are used to induce 
numerical damping and thereby reduce oscillation of the solution in 
regions of high flow compression. The method of characteristics first 
transforms the system of governing equations into characteristic form, 
after which the derivatives in the resulting equations are replaced by 
finite differences. The system of difference equations is then solved 
to obtain the solution. The advantages and disadvantages of each of 
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these methods have been discussed by Strom (1), Sauerwein (2), 

Rlchtmyer and Morton (3), and Forsythe and Wasow (4). Summarizing 
their findings, the finite difference methods are conceptually simpler, 
less difficult to program, require less computer storage, and can 
obtain the solution on an evenly spaced grid. The characteristics 
methods are, generally speaking, more accurate due to their more 
rigorous treatment of the physics of the problem. 

In the present investigation, the flow field is computed using 
the method of characteristics for steady three-dimensional flow. The 
bow shock wave and the internal shock wave system are computed using a 
three-dimensional discrete shock wave fitting procedure. Moreover, 
the influence of molecular transport may be included in the computation 
by treating the viscous and thermal diffusion terms in the governing 
equations of motion as forcing functions, or correction terms, in the 
method of characteristics scheme. The primary purpose in including the 
effects of molecular transport in the computation is for the possible 
future matching of the present analysis with a higher-order boundary 
layer analysis. No attempt was made in the present investigation to 
compute the boundary layer, or to account for boundary layer removal. 

3. GENERAL FEATURES OF THE THREE-DIMENSIONAL METHOD OF CHARACTERISTICS 

Extensive literature surveys of the method of characteristics 
for three-dimensional flow have been given by Zucrow and Huffman (5), 
Fowell (6), Thompson (7), Chushkin (8), Strom (1), Sauerwein (2), and 
Ransom, Hoffman, and Thompson (9). A brief summary of their conclu- 
sions is given here. 
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In genera], characteristics schemes for steady three-dimensional 
flow may be classified as either reference plane methods or bi char- 
acter! stic methods. In reference plane methods, the system of 
governing partial differentia] equations in three-independent variables 
is reduced to a system of partial differential equations in two inde- 
pendent variables by suitably approximating the derivatives with re- 
spect to the third independent variable. These approximations to the 
derivatives are then treated as forcing terms, and the resulting 
system of equations is solved using a two-dimensional characteristics 
scheme. Reference plane methods have been proposed by Ferrari (10), 
Sauer (11,12), Ferri (13), Moretti, et al. (14,15), Katskova and 
Chushkin (16), Holt (17,18), and Rakich (19). Reference plane methods 
have been called the method of bycharacteris vies by Moretti, et al. 
(14,15), the method of near characteristics by Sauer (11), the method 
of secondary characteristics by Sauer (12), and the method of semi- 
characteristics by Chushkin (8). In bicharacteristic methods, the 
characteristic equations are solved along the actual generators (bi- 
characteristics) of the Mach conoid and along the streamlines. Bichar- 
acteristic schemes have been proposed by Thornhill (20), Fowell (6), 
Sauerwein (2,21), Coburn and Dolph (22), Holt (23), Strom (1), Butler 
(24), and Cline and Hoffman (25). 

Reference plane methods are conceptually the simpler of the two 
schemes. However, reference plane methods have questionable accuracy 
in highly three-dimensional flows s‘>nce the domain of dependence of 
the differential equations is not rigorously considered. Al ternatively, 
while the bicharacteristic methods more rigorously treat the domain 
of dependence, they are also more complicated. The bi character! stic 
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methods are potentially the more accurate and, therefore, a bi char- 
acteristic method was selected for use in the present investigation. 

The particular bi characteristic method selected was that devised by 
D.S. Butler (24). Butler's scheme has been applied by Elliott (26), 
Richardson (27), and Delaney (28) to compute unsteady two-dimensional 
flows. Ransom, Hoffman, and Thompson (9) applied Butler’s method to 
compute the continuous steady three-dimensional supersonic isentropie 
flow field in nozzles, and Cline and Hoffman (25) applied Butler's 
method to compute the continuous steady three-dimensional supersonic 
flow field in nozzles accounting for nonequilibrium chemical reactions. 

A detailed description of the computer program developed for this 
calculation is given in NASA TM- 78947, "A Computer Program for the 
Calculation of the Flow Field in Supersonic Mixed- Comp res si on Inlets 
at Angle of Attack Using the Three-Dimensional Method of Characteristics 
with Discrete Shock Wave Fitting" by Joseph Vadyak, Joe D. Hoffman, and 
Allan R. Bishop. 
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SECTION II 
GAS DYNAMIC MODEL 


1. INTRODUCTION 

The gas dynamic model is based on the following major assump- 
tions : 

1. continuum flow 

2. steady flow 

3. negligible body forces 

4. the working gas can be represented as a simple system in 
thermodynamic equilibrium 

5. no mass diffusion 

6. negligible radiative heat transfer and no internal heat 
generation other than viscous dissipation 

7. viscous and thermal diffusion effects of secondary importance 
The governing equations for the assumed flow model consist of the con- 
tinuity equation, the component momentum equations, the energy equa- 
tion, the thermal and caloric equations of state, and appropriate 
representations for the molecular transport properties. These equa- 
tions are briefly presented in this section. A detailed development 
of these equations is given in Appendix A. 
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2. governing differential equations of motion 

it 

The continuity equation [see Reference (29}] is given by 


Dp 

Dt 


+ 


3u. 



= 0 


0 ) 


where x^ (i=l,2,3) denotes the three rectangular cartesian coordinates 
x, y, and z, repsectively , (i=l,2,3) denotes the corresponding 

velocity components u, v, and w, respectively, p denotes the density, 
and t denotes the time. The operator 0( }/0t in equation (1) is the 
material derivative given by 


MJU iLU 

Dt 9t 


£ LA 

9x . 


( 2 ) 


For steady three-dimensional flow, equation (1) may be written in ex- 
panded form as 


pu x + pVy + pw z + up x + vPy + wp^ = 0 (3) 

where the subscripts x, y, and z denote partial differentiation with 
respect to the corresponding direction. 

The momentum equation is given by the Navier-Stokes equation (29), 


which written in component form is 
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Repeated indices imply summation over the range of 1 to 3 unless other- 
wise noted. 





where EL (1=1 ,2,3} denotes the x, y, and z components of the body force, 
respectively, P denotes the pressure, p denotes the dynamic viscosity, 
and n denotes the second coefficient of viscosity. 

One of the major assumptions of the present investigation is 
that the influence of molecular transport is considered to be of 
secondary importance as compared to the inertial effects in determining 
the solution. As a consequence, the viscous and thermal diffusion 
terms appearing in the governing partial differential equations are 
treated as forcing functions or correction terms in the method of 
characteristics scheme to be presented. In what follows, the molecular 
transport terms are placed on the right-hand sides of the respective 
governing equations, and the convective terms are placed on the 
left-hand sides of those equations. The convective terms then are 
considered as constituting the principal parts of these equations. 

Hence, by assuming steady flow, negligible body forces, n = 0 [Stokes's 
hypothesis (30)], inertial dominance, and variable transport properties, 
equation (4) may be written in expanded form for each of the three 
coordinate directions as 


puu x + pvu y + pwu z + p x - F x 


puv x + pvv y + pwv z + P y - F y 


puw x + pvw y + pww z + ? z = F z 


(5) 

( 6 ) 
(7) 


where 


F„ - p 

+ P 


- p id U - l(v + w ) 
x x[3 x 3 y z_ 


+ V U y + \> + ^< U Z + 


4 u + u + u + l(v + w ) 

3 xx yy zz 3 V xy xz' 


( 8 ) 
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l (u x + w z^J + p x (v x * u y ) + p z (v z + w y> 

+ "6 v yy * v *x + v zz + i(“yx + "yz j] 

F z * u z[j w z ' f( u x + v y>] + “xK + «*> + “y^y + v z> 

+ “If w z* + "w + w yy + I (u zx + (10) 

The appropriate form of the energy equation is now derived. In the 
following, the pressure P and density p are considered as being the 
primary thermodynamic variables. All secondary thermodynamic vari- 
ables are then expressed in terms of the pressure and density. 

It is assumed in the present investigation that the working gas 
may be represented as a simple system in thermodynamic equilibrium. For 
a simple system, specification of any two independent thermodynamic 
properties defines the thermodynamic state of the system (31). Hence, 
the following functional relationship may be written 

P = P(p,s) (11) 




v - 


where s is the entropy per unit mass. Employing the concept of the 
total derivative, and introducing the material derivative operator 
given by equation (2), the following equation is obtained. 


DP 

f|P| 

te + 

fifl 

Dt " 

M 

s Dt 
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The sonic speed a is defined by 


Ds 

Dt 
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( 12 ) 


(13) 
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Introducing equation (13) into equation (12) yields 


0E . a 2 & = f2£l Pi 

“ l 8s J Dt 


Dt 


Dt 


(14) 


The material derivative of entropy appearing in equation (14) may 
be expressed in terms of a thermal conduction function and a viscous 
dissipation function. The entropy may be expressed in terms of the 

4 

internal energy by use of the thermodynamic relation (31) 


T ds * de + P d(l/p) 


05) 


where T is the temperature, and e is the internal energy per unit mass. 
The internal energy may be expressed in terms of a thermal conduction 
function and a viscous dissipation function by use of the energy 
equation (29) 


De _ 3 L 3T 1 


p Dt = 3x 


i 


K ~T— ! +—!■?+$' 


ax i) 


P Dt 


(16) 


where k is the thermal conductivity, and $ is the viscous dissipation 
function, which for n = 0 is given by 

au/ 


. 1 
$ = 2 u 


-.2 


ax. 


ax 


2 3u k ■ 

3 ax k °ij 


(17) 


where 6^ is the Kronecker delta. Combining equations (14) to (17) and 
writing the resulting expression in expanded form for steady three- 
dimensional flow with variable transport properties yields 


UP X + vP y + wP z - a (up x + v Py + wp z ) = F e 


(18) 


where 
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= s{k(T+T + T J + k T + k T H T 
l ' xx yy zz ' :< x y y z z 
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y x 


Z X 


z y' 


+ W x + u y + w y + u z + *2 - l (u x + v y + < 19 > 



( 20 ) 


3. THERMODYNAMIC MODEL 

Before a solution to the system of governing partial differential 
equations may be obtained, the temperature T, the sonic speed a, and x 

the parameter 5 defined by equation (20) must be expressed in terms 
of the primary thermodynamic variables P and p. The general functional 
forms of the relations for T, a, and f; are given by 


T = T(P,p) 

(21) 

a = a(P,p) 

(22) 

€ - C(P.p) 

(23) 


The derivatives of the temperature appearing in equation (19) are ex- 
pressed in terms of the derivatives of the pressure and the density 
by analytically differentiating equation (21). 

For the special case of a thermally and calorically perfect gas, 
equations (21) to (23) take the following simple forms 


T = P/pR 

(24) 

a = (yP/p) 1/Z 

(25) 




K • Y - 1 


(26) 


where R is the gas constant, and y is the specific heat ratio. 
4. MOLECULAR TRANSPORT PROPERTIES 


The dynamic viscosity v and the thermal conductivity < must be 
expressed in terms of the primary thermodynamic variables P and p. In 
general, both the viscosity and the thermal conductivity are assumed 
to be functions of temperature only. Hence, 


y * y(T) 

K = tc('f) 


( 2 ?) 

(28) 


The derivatives of the transport properties appearing in equations (8), 
(9), (10), and (19) are obtained in terms of the derivatives of the 
pressure and the density by analytically differentiating equations (27) 
and (28) with respect to the temperature, with the resulting tempera- 
ture derivatives being obtained by analytically differentiating equa- 
tion (21). 

A widely accepted representation for equation (27) is the 
Sutherland formula (30) 


v* = K 



1.5 

c._ _ 

T + S 

T 


0 

T 


T + S 

oj 


L 


(29) 


where u 0 is the viscosity at the reference temperature T Q » and S is a 
constant. Equation (28) may be represented by the quadratic polynomial 

r 2 


ic = a.j + a^T + a^T 


(30) 
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where the coefficients a^ (i=l,2,3) are obtained by curve fitting thermal 
conductivity data. 

The contribution of turbulent transport may be considered in the 
computation by adding the appropriate eddy viscosity and eddy thermal 
conductivity functions to the molecular transport properties given by 
equations (27) and (28), respectively. 

5. SUMMARY 

In summary, the differential equations of motion for steady three- 
dimensional flow are given by equations (3), (5), (6), (7), and (18). 

For a thermally and calorically perfect gas, the thermodynamic model 
is represented by equations (24) to (26). The molecular transport 
properties are represented by equations (29) and (30). 
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SECTION III 


CHARACTERISTIC EQUATIONS 


1. INTRODUCTION 

Written in the form shown, with the left-hand sides constituting 
the principal parts, equations (3), (5), (6), (7), and (18) may be 
classified as a system of quasl-linear nonhomogenous partial differ- 
ential equations of first order. The system is hyperbolic if the 
flow is supersonic. Systems of hyperbolic partial differential equa- 
tions in three independent variables have the property that there 
exist surfaces in three-dimensional space on which linear combinations 
of the original partial differential equations can be formed that con- 
tain derivatives only in the surfaces themselves. These special sur- 
faces are known as characteristic surfaces, and the linear combinations 
of the original partial differential equations are interior differen- 
tial operators known as compatibility relations. In this section, the 
equations for the characteristic surfaces and the compatibility rela- 
tions valid along these surfaces are listed and briefly discussed. A 
detailed development of these equations is given in Appendix B. 

2. CHARACTERISTIC SURFACES 

For steady three-dimensional supersonic flow, two families of 
characteristic surfaces exist, as illustrated in Figure 2. One family 
of characteristic surfaces consists of the stream surfaces given by 
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FIGURE 2. CHARACTERISTIC SURFACES 




(31) ' 


UN + VN + wN * 0 
x y z 

where N = (N ,M ,N ) denotes the normal to a stream surface. The 
x y 2 

envelope of all stream surfaces at a point forms a single pencil of 
planes whose axis is a streamline. A streamline may be represented by 

dx/dt = u dy/dt = v dz/dt = w (32) 


where t is the time of travel of a fluid particle along the streamline. 

The second family of characteristic surfaces consists of the wave 
surfaces given by 


uN + vN + wN = a|N { (33) 

x y z 

where N = (N »N ,N ) denotes the normal to a wave surface. The envelope 
x y z 

of all wave surfaces at a point forms a conoid known as the Mach conoid. 
The Mach conoid may be represented locally by a right circular cone 
known as the Mach cone. In differential form, the quadric surface of 
the Mach conoid is given by 


[ii 2 - (q 2 - a 2 )](dx) 2 ♦ [v 2 - <q 2 - a 2 )J(dy) 2 
+ [w 2 - (q 2 - a 2 )](dz) 2 + 2uv(dx)(dy) 

+ 2uw(dx)(dz) + 2vw(dy)(dz) = 0 (34) 


2 2 2 2 

where q is the velocity magnitude (q - u + v + w ). The line of 
contact between a particular wave surface and the Mach conoid is known 
as a bicharacteristic. A bieharaeteristic is a generator of the Mach 
conoid. 


19 


■SgSSR ■ ■■taw 


-■ i 


Ktl3w¥^V.^.*!iltar4.*i^- rr r jr, „ — »i n t , .. , 




3. COMPATIBILITY RELATIONS 

The compatibility relations which are applicable on the stream 
surfaces are given by 

uP x + vP y + wP 2 * a 2 (uP x + vp y + wp z ) = F fi (35) 

pu{uu + vu + wu ) + pv(uv + vv + wv ) 

« j *> a y z 

+ pw(uw + vw + ww ) + uP + vP + wP 

A Jr 4 A Jr 4 

- uF x ♦ *F y + wF z (36) 

PS X ( U U X + vu y + wu z ) + pS y (uv x + vv y + wv z ) 

+ pS z (, % + + "%> + Vx + S y P y + S z P z 

■ S x F x * S y F y + Vz < 37 > 

In equation (37), S = (S »S ,S ) denotes a vector which lies in the 

x y z 

stream surface and that is independent of the velocity vector. Equa- 
tions (35) and (36) may be written in a form which contains differen- 
tiation in the streamline direction as follows. 

ft ' ft ■ F e (38> 

pu It + pv M * pw 3? + 3F * uF x + ^y * wF z (39) 
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In equations (38) and (39), the operator d( )/dt represents the direc- 
tional derivative along a streamline. 
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The compatibility relation which is applicable on the wave sur- 
faces is given by 


pan x (uu x + vu y + wu z ) + pan y (uv x + vv y + wv z ) 

f pan (uw + vw + ww ) + (an - u)P + (an - v)P 
z * y *■ * * y j 

+ (an z - w)P z - pa 2 (u x + v^ + w 7 ) = A (40) 

where 


A = a(n x F x + n y F y + n z F z ) - F fi (41) 

In equations (40) and (41), n = (n ,n ,n ) denotes the unit normal 

x y z 

vector to the wave surface. Equation (40) may be written in a form 
which contains differentiation in the bicharacteristic direction as 
follows. 


pan x Ilk pan y M * pan z " ill = A ' pa?t < n x ‘ 


? 7 

+ (n - 1 )v + (n - 1 )w + n n (u +■ v ) 
\ y y v z z x y' y x / 


+ n *" 2 (u Z + "> ) * "y n z (v 2 + w y )] 


(42) 


In equation (42), the operator d( )/dt denotes the directional derivati 
along a bicharacteristic. The terms in brackets in equation (42) repre 
sent differentiation in the wave surface but in a direction normal to 
the: bicharacteristic direction. Hereafter, these terms will be refer- 
red to as the cross derivatives. 

At any point in the flow field there exists an infinite number of 
stream surfaces and wave surfaces. The number of independent 
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compatibility relations cannot exceed the number of independent equa- 
tions of motion. As a consequence, it is necessary to determine which 
of the possible combinations of the compatibility relations form an 
independent set. Rusanov (32), using a proof in the space of char- 
acteristic normals, has shown for steady three-dimensional isentropie 
flow that two of the stream surface compability relations applied 
along a stream surface and the single wave surface compatibility 
relation applied along three different wave surfaces form an indepen- 
dent set of five characteristic relations. Rusanov's results may be 
extended to the present case since the principal parts of equations 
(3), (5), (6), (7), and (18) are the same as those for isentropie 
flow. Hence, the set of compatibility relations used in the present 
investigation consists of equations (38) and (39) applied along a 
streamline and equation (42) applied along three different bi character- 
istics. 


4. BUTLER’S PARAMETERIZATION OF THE CHARACTERISTIC EQUATIONS 

D. S. Butler (24) developed a parameteric form for representing 
a bi characteristic and the wave surface compatibility relation 
applicable along it. A detailed development of Butler’s method is 
presented in Appendix B. A brief summary is given here. 

Butler introduced the following parameteric form to represent a 
bi characteristic. 

dx. = (u i + ca.cose + c^sin6)dt (i=l,2,3) (43) 

In equation (43), t is the time of travel of a fluid particle along 
the streamline that is the axis of the Mach cone, P is a parametric 
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angle denoting a particular element of the Mach cone and has the range 
0 <_ 0 <_ 2tt, and c is given by 

c 2 = q 2 a 2 /(q 2 - a 2 ) (44) 

where q is the velocity magnitude, and a is the sonic speed. The 
vectors a. and in equation (43) are parametric unit vectors with 
B^» and u^/q (1=1 ,2,3) forming an orthonormal set. 

The corresponding parametric form of the wave surface compatibility 
relation, equation (40), is given by 

HP du i 

g- + pc{«j cos6 + Basins) = $ 

? 3u i 

- pc (cx.sin© - g.cos0)(a,sin© - 3 .cose) t— (45) 

i l j j ax j 

In equation (45), the operator d( )/dt represents differentiation in 
the bi characteristic direction, and $ is given by 

4 = (c 2 /a 2 )[F e - a(n x F x + n y F y + n^)] ( 46 ) 

A 

where n = (n ,n ,n ) denotes the unit normal to the wave surface, 
x y z 

which may be written in parametric form as 

n,. = (a/c)(cu../q 2 - a. cose - B^ine) (i=l,2,3) (47) 

In addition to the above relations, Butler also developed a non- 
characteristic relation which is applied along a streamline. This 
noncharacteristic relation is given by 


(48) 


dP 

dt 


= a - 


P c ( a j^j 


+ $.$.) 

P 1 T 



where the operator d( )/dt denotes differentiation along a streamline, 
and (j is given by 

a - (c 2 /a 2 )F e - (c 2 /q 2 )(uF x + vF y + w? 2 ) (49) 


* 
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SECTION IV 


UNIT PROCESSES 


1. INTRODUCTION 

A variety of unit processes are employed in the computation of 
the flow field. The unit processes may be classified into four major 
types: interior point, solid boundary point, field-shock wave point, 

and solid body-shock wave point. The basic unit processes are briefly 
discussed in this section. A detailed presentation of each unit 
process is given in Appendix E. 

In the overall numerical algorithm, an inverse marching scheme is 
employed. The solution is obtained on space-like planes of constant x, 
where the x-axis is the longitudinal axis of the centerbody and the 
cowl. For the internal flow, the solution is also obtained on the 
space curves which are defined by the intersections of the internal 
shock wave with the solid boundaries. Except in the vicinity of a 
shock wave- sol id boundary intersection, the distance Ax between suc- 
cessive solution planes is determined by the application of the 
Courant-Friedrichs-Lewy (CFL) stability criterion (9). In the vicinity 
of a shock wave-solid boundary intersection, the axial step is con- 
trolled by special constraints, which are discussed in Section V. The 
distance Ax is determined prior to the application of the unit pro- 
cesses. 
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2 . INTERIOR POINT UNIT PROCESS 

The computational network used in determining the solution for a 
typical interior point is illustrated in Figure 3. Points (1) to (4) 
represent the intersection points of four rearward-running bicharacter- 
isties with the initial- value plane, point (5) is the streamline 
intersection point with the initial-value plane, and point (6) is the 
solution point on the solution plane. The axial (x) distance between 
the initial-value plane and the solution plane is determined prior to 
the application of the unit process by applying the CFL stability 
criterion. As in all the unit processes, the interior point unit pro- 
cess is divided into a predictor step and a corrector step. The cor- 
rector may be iterated to convergence if desired. 

The interior point unit process is initiated by determining the 
location of the solution point, point (6). The coordinates of point 
(6) are determined by extending the streamline from point (5) to the 
solution plane using the following finite difference form of equation 
(32). 

x. (6) - Xj (5) = ]- [ Uj (5) + Uj (6)1[t(6) - t{5)3 (1-1.2.3) (50) 

For the predictor, u^(6) is equated to u. (5). For the corrector, 
the previously determined value of u^(6) is employed. The axial step 
[x{6) - x(5)] is computed before the unit process is applied. Hence, 
the time parameter [t(6) - t (5 ) ] may be obtained, after which the 
coordinates y (6) and z(6) are computed. Interpolated flow property 
values at point (5) are used in the integration, even though point (5) 
is a known field point. As shown by Ransom, Hoffman, and Thompson (9), 
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FIGURE 3. INTERIOR POINT COMPUTATIONAL NETWORK 


no 


this interpolation is required to produce a stable numerical scheme. 

The interpolated flow property values are obtained from the following 
quadratic bivariate interpolation polynomial 

2 2 

f(y,z) ■ + agy + a 3 z + a 4 yz + a g y + a g z (51) 

where f(y,z) denotes a general function of the coordinates y and z, and 
the coefficients a^ ( i =1 to 6) are obtained from a least squares fit 
of nine data points in the initial -value plane [point (5) and its 
eight immediate neighbors] as described in Appendix C. 

With the location of the solution point determined, four bi char- 
acteristics are extended from the solution point back to the initial- 
value plane to intersect this plane at points (1) to (4), as illustrated 
in Figure 3. The coordinates of each of these intersection points are 
determined using the following finite difference form of equation (43). 

x. (6) - (k ) = j {u. (k) + u^ (6) + [c(k) + c(6)][a.cos6(k) 

+ 0jSln6(k)]}[t(6) - t(k)] (i=l ,2,3) (52) 

The index k in equation (52) denotes the bi characteristic- initial - 
value plane intersection points illustrated in Figure 3, and has a 
range of 1 to 4, corresponding to the 9(k) values of 0, ir/2, tt, and 
3 tt/ 2, respectively. Since the axial step [x(6) - x(k)] is known, 
equation (52) is used to calculate [t(6) - t(k)], y(k), and z(k). The 
flow property values at points (1) to (4) are obtained by interpolation 
using equation (51). On the initial application of equation (52), the 
flow property values at point (k) are equated to those at point (5). 
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For the external flow field integration, the parametric unit vectors 
a. and 8. appearing in equation (52) are selected to straddle the pro- 
jection of the pressure gradient on the initial-value plane. For the 
internal flow field integration, these vectors are selected to straddle 
the meridional plane through point (6). 

Once the positions of and the flow properties at points (1) to 
(S) have been determined, the system of nonlinear compatibility equa- 
tions, written in finite difference form, is solved to obtain the five 
dependent flow properties u(6), v(6), w(6), P {6 ) and p(6). Two of the 
five required compatibility equations are given by equations (38} and 
(39). These equations are written in finite difference form by re- 
placing the derivatives with simple differences, and by replacing the 
coefficients of the derivatives with the arithmetic average of the 
coefficients at the solution point and at the appropriate point in the 
initial-value plane. To obtain the remaining three required compati- 
bility equations, appropriate linear combinations of the wave surface 
compatibility relation, equation (45), applied along each of the four 
bi characteristics, and the noncharacteristic relation, equation (48), 
applied along the streamline are formed. Writing equation (45) for 9 
values of 0, it/ 2, tt, and 37t/2 yields 


dP 


du.. 


dt 1 + pCGt i dt 1 * *1 " pC 6 i^j 


3u^ 

ax' 


(53) 


dP_ + 
dt 2 


du. 


1 . * 


9U. 

1 


i * *2 - pc a i“j 3Xj 


(54) 




In equations (53) to (56), the operator d( )/dt^ denotes differentiation 
along the kth ^characteristic, and denotes equation (46) evaluated 
for the specified value of © (k ) . One independent linear combination 
of the compatibility equations is obtained by subtracting the finite 
difference form of equation (55) from the finite difference form of 
equation (53). Another independent linear combination is obtained by 
subtracting the finite difference form of equation (56) from the finite 
form of equation (54). The final independent linear combination is 
obtained by subtracting the finite difference form of the nonchar- 
acteristic relation, equation (48), from the sum of the finite differ- 
ence forms of equations (53) and (54). The resulting compatibility 
equations do not contain cross derivatives at the solution point [i.e., 
all terms containing 3u./3xj(6) are eliminated]. These five finite 
difference equations are solved using Gaussian elimination. For the 
predictor, the flow property values at the solution point appearing 
in the coefficients of the derivatives in the set of difference equa- 
tions are equated to those at point (5). For the corrector, the flow 
property values at point (6) obtained on the previous iteration are 
used. The resulting scheme has second-order accuracy (9). 
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3. SOLID BOUNDARY POINT UNIT PROCESS 

The computational network used for determining the solution at a 
typical point on a solid boundary is shown in Figure 4. The point 
notation used in this figure is identical to that employed in Figure 3. 
Here, however, both points (5) and (6) lie on the solid boundary, and 
point (4) is not used since it lies outside of the flow regime. 

The unit process used to obtain the solution at a solid boundary 
point is almost identical to the interior point unit process. Here, 
however, point (4) corresponding to the bicharacteristic with e = 3tt/2 
is not located, and the corresponding compatibility relation valid 
along this bi characteristic is not employed. That equation is replaced 
by the flow tangency condition 

u i ( 4 * 6 ) n bi (6) = 0 (57) 

where n bl . (6 ) (i=l,2,3) is the unit normal to the solid boundary at 
point (6). 

4. BOW SHOCK WAVE POINT UNIT PROCESS 

The computational network used in determining the solution for a 
typical bow shock wave point is illustrated in Figure 5. A segment 
of the shock wave surface extending from the initial-value plane to 

the solution plane is shown in this figure. The intersection of the 
shock wave with the initial-value plane defines space curve (A), and 
the intersection of the shock wave with the solution plane defines 
space curve (B). The axial distance between the initial-value plane 
and the solution plane has been previously determined by the applica- 
tion of the CFL stability criterion. The bow shock wave solution point 
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FIGURE A. SOLID BOUNDARY POINT COMPUTATIONAL NETWORK 




WAVE POINT COMPUTATIONAL NETWORK 



is denoted by point (2). The flow properties at point (2) on the up- 
stream side of the shock wave are known from the free-stream condi- 
tions. Hence, in the following discussion, the flow properties u(2), 
v(2), w(2), P(2), and p{2) refer to the flow properties at point (2) 
on the downstream side of the shock wave. Point {1) is the intersec- 
tion point of a rearward -running bi characteristic with the initial- 
value plane. This bi characteristic is extended backward from the 
solution point, point (2). Point (3) is a predetermined interior solu- 
tion point which is adjacent to the shock wave and is used to define 
the meridional plane in which the bow shock wave solution point lies. 
Point (4) is the intersection point of space curve (A) with the 
meridional plane which passes through points (2) and (3). 

In this unit process, a local cartesian coordinate system is 
employed for the description of the local shock wave surface. This 
local coordinate system has cocrdinates x', y‘, and z', where x' is 
coincident with the x-axis, y* is the radial direction in the meridional 
plane containing points (2) and (3), and z’ is normal to the (x'.y')- 
plane. The unit vectors in the x', y’, and z' directions are denoted 

A A A 

by i‘, j'» and k*, respectively. The orientation of the local shock 
wave surface at a point (P) is specified by a set of three unit vec- 
tors referenced to the (x' »y 1 ,z' ^coordinate system, as illustrated 

A 

in Figure 6. This set of unit vectors consists of the unit vector n s 
which is normal to the shock wave surface at point (P), and two unit 

A A 

vectors £ and t which are tangent to this surface at point (P). The 

A 

tangential unit vector t lies in the meridional plane [{x* ,y’ )-plane], 
subtends an angle $ with the x'-axis, and is defined by the intersection 
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of the shock wave with the meridional plane at point (P). The tan- 

a 

genual unit vector £ lies in the transverse plane [ (y * ,z' )-plane], 
subtends an angle n with the z'-axis, and is defined by the intersec- 
tion of the shock wave with the transverse plane at point { P ) . The 

a. 

tangential unit vectors t and t are given by 

a A A 

t = cos i 1 + sin $ j ' (58) 

J\ ■■A 

l ~ sin a j ' + cos a k“ (59) 

A 

The shock wave normal unit vector n s is given by 

n = H x t/jt t| (60) 

To achieve second-order accuracy in the shock wave point unit 
process, global iteration must be performed. In global iteration, 
the corrector employs flow properties not only at the solution point 
itself, but also at neighboring points in the solution plane. As a 
consequence, before the corrector can be applied in global iteration, 
the entire solution plane (or at least an appropriate section of it) 
must be determined by a prior calculation. The interior point and 
solid boundary point unit processes do not require global iteration 
to achieve second-order accuracy. Consequently, those solution points 
are determined first. Then, the predictor is applied for each shock 
wave solution point, thereby giving a tentative solution for all of 
the shock wave points. At this stage, global correction is performed 
for the shock wave solution points using the previously determined 
field points in the solution plane. In the following discussion, the 
term "predictor" refers to the first application of the shock wave 
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point unit process used to obtain an initial estimate of the solution 
without using field point data in the solution plane. The term "global 
corrector" refers to the application of the shock wave point unit 
process which uses field point data in the solution plane. The shock 
wave point unit process is now outlined. 

The shock wave point unit process is initiated by locating the 
solution point, point ( 2 ) in Figure 5. Denote the angle subtended by 
a meridional plane and the (x,y)-plane by 0. The solution point 
meridional plane is arbitrarily selected to contain the interior solu- 
tion point, point (3), whose location is determined prior to the 
application of the shock wave point unit process. Hence, ©(2) - 6(3). 
Denote the radial position of a point by r. Then the radial position 
of point (2) is obtained from 

r(2) = r(4) + [x(2) - x(4)] tan |^[<J>{2) + <f>(4)]j (61) 

where [x{2) - x(4)] is the axial distance between the initial-value 
plane and the solution plane. On the initial application of equation 
(61), the shock wave angle (f)(2) is equated to $(4), whereas, on ensuing 
applications, the value of (f>(2) obtained on the previous iteration is 
used. At point (4), the radial position r(4) and the shock wave angle 
(f)(4) are determined by interpolation using the quadratic univariate 
formulae 


r(e) 

2 

= aj + a 2 © + a 3 © 

(62) 

(f)(0) 

= b-j + b 2 © + b 3 © 2 

(63) 
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where the coefficients a^ and b.. { i -7 *2,3) are determined by fitting 
these expressions to three local shock wave solution points on space 
curve (A), 

After the solution point has been located, the shock wave normal 

A 

unit vector n $ at the solution point is found by forming the normal- 

A /S 

ized cross product of the tangential unit vectors I and t [see equation 

A J\ 

(60)]. The tangential unit vectors t and i are obtained by use of 
the current values of <f>(2) and a (2) in equations (S8) and {59), 
respectively. For a predictor application, a(2) is approximated by 
equating it to the a value at point (4). For a global corrector 
application, the value of ot{2) that is employed is that evaluated at 
point (2). In either case, the value of a(2) may be determined by 

«(2) = tan- T (i|^| 8(2) (64) 

where, for the predictor, the analytical form of r(e) used in equation 
(64) is given by equation (62) applied along space curve (A), and for 
the global corrector, r(e) is obtained by applying equation (62) along 
space curve (B). 

At this stage, the local Hugoniot relations are applied at point 
(2) to obtain the downstream flow properties u(2), v(2), w(2), P{2), 
and p{2). Next, a rearward-running bicharacteristic is extended from 
the solution point, point (2), back to the initial-value plane, inter- 
secting this plane at point (1), as illustrated in Figure 5. The 
coordinates of point (1) are obtained using the following finite 
difference form of equation (43) evaluated for the parametric angle of 
9 = %/Z. 
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( 65 ) 


(2) - Xl (1) = | <u 1 (1 } + U 1 (2) 

+ [c(l) + c(2)]B.) [t(2) - t(l)] (i=l ,2,3) 

For the first application of equation (65), the flow properties at 
point (1) are equated to those at point (2), whereas, for ensuing 
applications, the flow properties previously obtained at point (1) are 
employed. The flow properties at point (1) are obtained by interpola- 
tion using the quadratic bivariate polynomial given by equation (51). 
Since the axial step [x(2) - x(l)] is determined by the CFL stability 
criterion, equation (65) is used to compute [t(2) - t(l )], y (1 ) , and 
z{l). The orientation of the parametric vector in equation (65) is 
selected so that this vector lies in the meridional plane that con- 
tains the solution point. The unit vector au is obtained using the 
orthonormal relationship between ou, , and u^/q (i=l,2,3). 

At this stage, the wave surface compatibility equation correspond- 
ing to the parametric angle © - tt/ 2 is applied between points (1) and 
(2). The appropriate equation is obtained by writing equation (54) in 

finite difference form and solving for the pressure at point (2). De- 

★ 

note this pressure by P (2). The resulting equation contains cross 
derivatives (terms containing 3u./8x.) at both points (1) and (2). 

* J 

For the predictor, the cross derivatives at point (2) are equated to 
those at point (1), whereas, for the global corrector, the cross 
derivatives at point (2) are evaluated at that point by fitting inter- 
polation polynomials in the solution plane. 

The pressure P(2) is calculated from the local Hugoniot equations. 
The pressure P*(2) is calculated from the wave surface compatibility 
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relation. The difference between P{2) and P (2) is driven to within 
a speci fifed tolerance of zero using the secant method with iteration 
being performed on the shock wave angle <t>{2). Two initial estimates 
of f(2) are required to start the iterative process. 

The shock wave point unit process is first applied as a predictor 
for each shock wave solution point. In this application, the value of ^ 

a used in equation (59) is obtained by curve fitting points along space 
curve (A), and the cross derivatives at the solution point are equated 
to those at the bicharacteristic base point in the initial-value plane. 

After a tentative solution has been obtained at each shock wave point, 
a number of ensuing global corrections are performed. Here, the value 
of a used in equation (59) is based on data along space curve (B), and 
the cross derivative terms at the solution point are evaluated at that 
point. The resulting overall algorithm has second-order accuracy when 
the global correction is performed. The global iteration is terminated 
when successive values of a have converged at each of the shock wave 
solution points. 

5. SOLID BODY-SHOCK WAVE POINT UNIT PROCESS 

The solid body-shock wave point unit process is used to determine 
the flow properties downstream of the shock wave at a point where the 
shock wave intersects a solid boundary. This unit process is used to 
determine the solution for the points on the cowl on the downstream 
side of the cowl lip shock wave, and for the points on the centerbody 
or cowl on the downstream side of an internal reflected shock wave. 

The method of computation is essentially the same for either applica- 
tion, For the internal shock wave reflection, the flow properties 


40 


downstream of the incident shock wave, which constitute the upstream 
flow properties for the reflected shock wave, are computed by the 
modified field-shock wave point unit process discussed in Appendix E. 

A depiction of the computational network used in the solid body- 
shock wave point unit process is presented in Figure 7. A typical 
solid body-shock wave solution point is denoted by point (P), with the 
outward unit normal vector to the solid boundary at this point denoted 

a 

by n^. The locus of solid body-shock wave solution points represents 
the intersection of the shock wave with the solid boundary and defines 
space curve (A) in Figure 7, The intersection of the shock wave with 
the meridional plane passing through point (P ) defines space curve (B), 
The unit vectors tangent to space curves (A) and (B) at point (P) are 

A A 

denoted by £ and t, respectively. The unit vector normal to the shock 

-A 

wave at point (P) is denoted by r, s . 

A 

As for the bow shock wave point unit process, the unit vectors £, 

A A 

t, and n are referenced to the local coordinate system (x'.y'.z'), 
where x', y', and z' have the same definitions as noted before. More- 

A 

over, the tangential unit vector t again lies in the meridional plane and 
is defined by equation (58). In this scheme, however, the tangential 

A 

unit vector £ does not lie in the (y ' ,z‘ )-plane in most cases, but 
rather can have a nonzero x' -component. This tangential unit vector 
along space curve (A) may be represented by 


£ 


= dx* f . + *L' 
ds 1 ds 



I 


where ds is the differential arc length given by 
(ds) 2 = (dx 1 ) 2 + (dy') 2 * (dz') 2 


( 66 ) 


(67) 
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SPACE CURVE (A) 

FIGURE 7. SOLID BODY- SHOCK WAVE POINT COMPUTATIONAL 


NETWORK 



The derivatives in equation (66) are obtained by analytically differ- 
entiating the expressions 


x'(0) = a 1 + 

a 2 © 

♦ a 3» 

(63) 

y ' (© ) * b 1 + 

b 2 e 

+ <>/ 

(69) 

Z'(B) = C T + 

c 2 8 

+ c 3 e 2 

(70) 

where coefficients a.. 

V 

and c. 

(1=1 ,2,3) are obtained by curve 


fitting the respective expressions to three points on space curve (A). 
For the cowl lip shock wave points, space curve (A) is defined by the 
cowl lip itself, since the shock wave is assumed to be attached to the 
cowl lip. Alternatively, for computing the downstream flow properties 
at a reflected internal shock wave, space curve (A) is defined by 
the intersection of the incident shock wave with the solid boundary. 

The shock wave normal unit vector is found from equation (60), 

The solid body-shock wave point unit process is initiated by 
determining the body normal unit vector n^ and the tangential unit 
vector JL An assumption is then made for the shock wave angle <p In 
equation (58), and, by use of equation (60), the shock wave normal unit 
vector is determined. The local Hugoniot equations are then applied to 
obtain the downstream flow properties at point (P). The velocity 
normal to the wall is then obtained by forming the dot product of the 
body normal vector and the downstream velocity vector. The normal 
velocity is reduced to within a specified tolerance of zero by varying 
the shock angle $ using the secant iteration method. 
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6. INTERNAL FLOW SHOCK WAVE POINT UNIT PROCESSES 

The unit process employed to compute the solution at a shock wave 
point in the internal flow field is similar to the bow shock wave point 
unit process. In the internal flow shock wave point unit process, 
however, the flow properties upstream of the shock wave at the solution 
point must be determined by the application of a modified interior 
point unit process. Moreover, modifications to the internal flow 
shock wave point unit process must be made when an internal flow shock 
wave solution point lies on or close to a solid boundary. The various 
versions of the internal flow shock wave point unit process are pre- 
sented in Appendix E. 

7. INTERNAL SHOCK MOD I FI ED- INTER IOR POINT AND -SOLID BODY POINT UNIT 
PROCESSES 

In some situations during the computation of the interna? flow 
field, the interior point and solid boundary point unit processes 
must be applied in a modified form. One such instance in which 
a modified form of the interior point unit process must be applied is 
shown in Figure 8. Here, the Mach cone, with apex at the interior 
solution point, intersects not only the initial-value plane but also 
the internal shock wave and a solid boundary. The unit process used in 
this case requires determining the bicharacteristic intersection points 
with the shock wave and the solid boundary in addition to the inter- 
section points with the initial-value plane. Moreover, flow property 
values must be determined at all of these points. The bieharacteristic- 
shock wave and bi characteristic-body intersection coordinates are cal- 
culated using the procedures discussed in Appendix D. The flow 



FIGURE 8- SHOCK- MODIFIED INTERIOR POINT COMPUTATIONAL 
NETWORK (STREAMLINE BASE POINT ON 
INITIAL- VALUE PLANE) 
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property values at these points are obtained by interpolation, either 
using a quadratic bivariate polynomial [equation (51)] for points on 
the initial-value plane, or using a quadratic trivariate polynomial 
for points on the shock wave surface or solid boundary surface. The 
various interpolation schemes are discussed in Appendix C. All of the 
unit processes, including the schemes incorporating the necessary 
modifications to handle the internal shock wave, are presented in 

■K 

Appendix E. 
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SECTION V 


OVERALL NUMERICAL ALGORITHM 


1 . INTRODUCTION 

The overall numerical algorithm consists of the repetitive appli- 
cation of the various unit processes to generate the global solution 
for given boundary conditions and a specified set of initial data. 

The contours of the centerbody and the cowl, in addition to any 
planes of flow symmetry, constitute the boundaries of the computational 
flow regime. For the external flow field integration, the bow shock 
wave also represents a computational bound. 

The initial data are specified on a plane of constant x. The x- 
coordinate axis is the longitudinal axis of the centerbody and the cowl 
(see Figure 1). Moreover, the mean flow direction is assumed to be in 
the x-coordinate direction. 

An inverse marching scheme is employed in the numerical algorithm. 
The solution is obtained on space- like planes of constant x. The solu- 
tion points on each plane represent the intersection points of continu- 
ous streamlines which are propagated from the data points specified on 
the initial-value plane. In addition to the streamline solution points, 
solution points are also obtained at the intersection of the external 
and internal shock waves with the solution plane, and for the internal 
flow field, on the space curves where the internal shock wave intersects 
the solid boundaries. These space curves are defined by the locus of 
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shock wave solution points. 

Except in the vicinity of a shock wave reflection with a solid 
boundary, the axial (x) distance between the current initial-value 
plane and the current solution plane is detained by the application 
of the Courant-Friedrichs-Lewy (CFL) stability, criterion (9). In the 
vicinity of a shock wave reflection with a soli a. boundary, the axial 
distance between successive solution planes is chosen so that the 
entire shock wave-solid boundary intersection falls between two 
adjacent solution planes. 

The external flow about the forebody is computed first. The ex- 
ternal flow field integration requires the periodic addition of 
streamlines in order to retain a well dispersed computational mesh. 
Furthermore, periodic deletion of selected streamlines is also re- 
quired so that the number of computational points lies within bounds. 

The internal flow field can be computed with or without the 
discrete fitting of the internal shock wave system. The option in which 
shock waves are not discretely fitted may be used in cases in which 
the internal shock waves are quite weak in strength, and thereby an 
acceptable solution can be obtained by smearing the internal discon- 
tinuities. 

In this section, brief discussions are presented on generation of 
the initial data, boundary conditions, regulation of the marching step 
size, computation of the transport forcing functions, and numerical 
stability. A detailed discussion of the overall numerical algorithm 
is presented in Appendix F. 
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2. INITIAL* VALUE PLANE 


The initial data are specified on a plane of constant x (see 
Figure 1). The flow must be supersonic at every point on this plane. 

For uniqueness and existence of a genuine solution* the values of the 
five dependent variables (u, v, w, P, and p) prescribed on this sur- 
face must have at least continuous first partial derivatives. 

If the forebody flow field is to be computed, the initial-value 
plane must be specified at an axial (x) station that is upstream of 
the forebody flow computational regime (see Figure 1). The last solu- 
tion plane of the forebody flow field computation is adjusted to lie at 
the axial station of the cowl lip, and constitutes the ini tial -value 
plane for the internal flow field computation. The cowl lip is assumed 
to be contained in a plane of constant x. Furthermore, the bow shock 
wave must fall outside of the cowl lip, or, in the limit, intersect 
the cowl lip exactly. The internal flow cannot be calculated if the 
bow shock wave is ingested into the annulus. The points on the solu- 
tion plane at the cowl lip axial station are redistributed to obtain a 
ring of solution points coincident with the cowl lip. 

If the forebody is conical ahead of the axial station where the 
initial-value plane is specified, an approximate flow property field 
on this plane may be internally generated in the computer program. The 
internally generated initial data are obtained by an approximate tech- 
nique which employs the Taylor-Maccol 1 solution for the flow about a 
circular cone at zero incidence. A superposition method is then used 
to obtain an approximation for the flow about a circular cone at 
nonzero angle of attack by neglecting the cross flow effects. Alterna- 
tively, a more exact solution for the initial data for flow about a 
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circular cone at incidence may be obtained by employing the results of 
Jones (33). 

If the forebody is not conical ahead of the axial station of the 
initial-value plane* another source of initial data must be used. If 
available, experimental data may be employed. 

3. SOLID BOUNDARY SURFACES 

The computer program developed in the present investigation assumes 
that both the eenterbody and the cowl are axi symmetric. For the pur- 
poses of geometry description, the axial (x) domain is divided into a 
number of intervals. In any interval, the body radius r may be speci- 
fied by either tabular input, or by supplying the coefficients in the 
cubic polynomial 

r(x) = + b. (x - ) + c..(x - ) 2 + d^(x - ) 3 (71) 

where the subscript i denotes the ith interval, r (x ) is the body radius 
at axial position x (x i <_ x < and the coefficients , b.. , 

and d. are obtained by curve fitting the body contour. Since equation 
(71) is a cubic, slope and curvature can be matched at the junction 
point between two adjacent intervals (i.e.» spline fits can be employed). 

4. INTEGRATION STEP SIZE REGULATION 

Except in the vicinity of a reflection of the internal shock wave 
with a solid boundary, the axial marching step between successive 
solution planes is determined by the application of the Gourant- 
Friedrichs-Lewy (CFL) stability criterion (9). The CFL stability cri- 
terion mandates that the domain of dependence of the differential equa- 
tions be contained within the convex hull of the finite difference 
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network. That is, the Mach cone must be inside the outer periphery of 
the nine initial-value plane field points used in formulating the bi- 
variate interpolation polynomial, equation (51). The allowable axial 
step is given by 

Ax = [u 2 /(cq)}[] - (c/q)(q 2 /u 2 - l) 1/2 ]R min < 7Z ) 

where Ax is the marching step, and R . is the distance between the 
streamline intersection point with the initial-value plane and the 
nearest point on the convex hull of the finite difference network. Equa- 
tion (72) is applied at every streamline point on the initial-value 
plane, with the actual integration step being chosen as the Ax 

value at the most restrictive point. Equation (72) is applied only to 
streamline points. The shock wave points are excluded. Moreover, in 
the internal flow field integration, the shock wave points are ignored 
in defining the convex hull of the finite difference network when 
applying the stability criterion to a streamline point. 

In the vicinity of a reflection of the internal shock wave with a 
solid boundary, the axial step is controlled by the constraint that the 
shock wave-solid body intersection is contained entirely between two 
adjacent solution planes. The fit point stencils used in formulating 
the various interpolation polynomials are appropriately expanded, in 
this case, so that the GFL stability criterion is satisfied. 

5. CALCULATION OF THE TRANSPORT FORCING FUNCTIONS 

The numerical procedure developed in the present investigation has 
the capability to include the influence of molecular transport on the 
solution by treating the viscous and thermal diffusion terms in the 
governing partial differential equations as forcing functions, or 
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correction terms, in the method of characteristics scheme. The computer pro- 
gram has the capability to include the influence of viscous and thermal diffusion 
in the computation of the external flow about the forebody, and in the computa- 
tion of the internal flow field in which shock waves are not discetely fitted. 
The program option in which shock waves are discretely fitted in the internal 
flow field does not have the capability to include the influence of molecular 
transport in the computation, but rather assumes the flow to be inviscid and 
adiabatic. The detailed calculation procedures used for obtaining the transport 
forcing terms are presented in Appendix G. 

6. NUMERICAL STABILITY 

A stability analysis of the nonlinear finite difference algorithm includ- 
ing molecular transport was not attempted. Instead, a stability analysis for 
isentropic flow was conducted. Stability of the generalized analysis was then 
verified by actual numerical calculations. 

Ransom, Hoffman, and Thompson (9) used the present numerical method to 
compute the continuous steady three-dimensional supersonic isentropic flow in 
a nozzle. The CFL stability criterion was used for locating successive solution 
planes, A von Neumann linear stability analysis indicated that interpolated 
flow properties, instead of the actual known values, should be used at the 
streamline-initial-value plane intersection point [point (5) in Figure 3]. The 
present analysis uses interpolated flow properties at all points in the initial- 
value plane. 


SECTION VI 


COMPUTATIONAL RESULTS 


1. INTRODUCTION 

Selected computational results are presented and discussed in this 
section. The results presented are divided into three major categories: 
external flow about the forebody, continuous internal flow, and in- 
ternal flow in which the internal shock wave system has been computed. 

In some instances, both axisynmetric flow and three-dimensional flow 
results are shown. For the internal flow field in which shock waves 
have been fitted, some comparisons with experimental data are made. 
Moreover, some results of the present analysis are compared with those 
of existing computational methods. 

2. EXTERNAL FLOW ABOUT THE FOREBODY 

For the purpose of testing the external flow integration procedure, 
the flow field about a right circular cone at incidence was computed. 

This is a conical flow field and the solution is constant along rays em- 
anating from the vertex of the cone (i.e., there is no characteristic 
length, so the solution has no dependency on x). At zero angle of 
attack, the solution depends only on the angle subtended by a given ray 
and the x-axis. At nonzero incidence, an azimuthal variation also 
exists. To obtain the required initial data, the results of Jones v.33 ) 
were employed. The computed results should maintain the conical nature 
of the flow field. 
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Figure 9 presents numerical results obtained for a 10.0° half- 
angle cone at 2.5° angle of attack a with a free-stream Mach number 
of 3.0. The computation employed HI circumferential stations in the 
computed sector (half-plane), and the number of radial stations on the 
initial -value plane was 11. The computed static pressure P normalized 
by the free-stream static pressure P m is plotted versus the axial 
position x normalized by the cowl lip radius R c . The pressure distri- 
butions on the rays formed by the forebody and the bow shock wave on 
both the leeward and windward planes of symmetry are shown. Since the 
flow is conical, the solution should remain constant along each of 
these four rays at the respective pressure values at the appropriate 
points on the initial-value plane. The initial-value plane pressures 
are denoted by the straight line segments. The method of characteris- 
tics solution is shown at a discrete number of axial stations, each 
station corresponding to the axial location of a given solution plane. 
The method of characteristics solution maintains the conical nature of 
the flow field. 

It should be noted that the increase in pressure across the leeward 
side of the bow shock wave is minimal. As the angle of incidence is 
further increased, the strength of the bow shock wave on the leeward 
side is reduced until the point is reached where the angle of attack 
is equal to the cone half-angle. At this point, no shock wave exists 
on the leeward meridional plane. Further increase in the angle of 
incidence causes a flow expansion to occur on the leeward side. Since 
the present analysis assumes that a shock wave exists about the 
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entire forebody, the case where a flow expansion occurs on the leeward 
side cannot be computed. 

Tie external flow about a circular cone at incidence was also 
computed including the effects of molecular diffusion. No significant 
changes in the computed results were noted. Approximately 60 percent 
more computer execution time was required for the computation which 
included the molecular diffusion terms. 

3. CONTINUOUS INTERNAL FLOW 

For the purpose of testing the continuous internal flow integra- 
tion procedure in which shock waves are not discretely fitted, the 
axisymmetric flow field in the simplified geometry inlet illustrated in 
Figure 10 was computed. The geometry of this inlet was selected so 
that the slope of the cowl contour at the cowl lip was equal to the 
slope of the streamline at the cowl lip. Hence, no flow turning occurs 
at the cowl lip and the internal shock wave system is not generated. 

Figure 10 illustrates the inlet geometry and the pressure dis- 
tributions on the centerbody and the cowl. A monotonic increase in 
pressure oh the surfa- e of the cowl occurs. The pressure on the center 
body retains its conical flow value until the Mach wave emanating from 
the cowl lip reaches the centerbody. After that point, a monotonic 
increase in the centerbody pressure occurs. This computation was per- 
formed with 21 radial stations and 1 circumferential station. The maxi 
mum deviation in mass flow rate on any solution plane from the mass 
flow rate across the cowl lip solution plane was 0.25 percent. 
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The flow field in the simplified geometry inlet illustrated in 
Figure 10 was also computed including the effects of molecular diffu- 
sion on the solution. No significant changes in the computed results 
were noted. The increase in computer execution time was approximately 
60 percent. 



4. INTERNAL FLOW WITH DISCRETE FITTING OF THE INTERNAL SHOCK WAVE 
SYSTEM 


Internal flow calculations were performed for the Boeing Mach 3.5 
supersonic mixed-compression inlet documented in Reference (34). The 
centerbody and cowl coordinates of this inlet are listed in Table 1. 

The boundary contours are illustrated in Figure 11 for the design case 
of zero centerbody translation. This inlet has a forebody which is 
conical (the forebody is not shown in Figure 11). Consequently, all 
of the numerical solutions were started at the cowl lip axial station. 
The initial data were obtained by employing the results of Jones (33). 

Extensive boundary layer removal is employed in this inlet to 
control boundary layer separation in regions of strong adverse pressure 
gradients such as those caused by oblique shock wave-boundary layer 
interactions. Figure 11 indicates regions where the boundary layer is 
removed. Since the present analysis does not compute the boundary la; "r 
nor takes account of boundary layer removal, good agreement between 
computed and experimental results cannot be expected in regions of high 
viscous interaction. For this inlet, 13.3 percent of the cowl lip mass 
flow rate was removed through boundary layer bleed at the design point 
condition to control boundary layer separation (34). 
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TABLE 1 


MACH 3.5 INLET COORDINATES 


CENTERBODY CENTER BODY COWL COWL 


x/R c 

r/R c 

x/R c 

r/R c 

x/R c 

r/R e 

x/R c 

r/R c 

0.0 

0.0 

4.8 

0.7504 

2.86 

1.0 

4.55 

0.8695 

4.0 

0.70532 

4.9 

0.7391 

3.1 

1.004188 

4.6 

0.864 

4.1 

0.7228 

5.1 

0.7120 

3.2 

1.0054 

4.65 

0.86 

4.2 

0.7387 

5.3 

0.6829 

3.4 

1.0051 

4.7 

0.8572 

4.3 

0.7512 

5.5 

0.6525 

3.6 

0.99996 

4.8 

0.8533 

4.4 

0.759 

5,6 

0.6362 

3.8 

0.9882 

4.9 

0.8511 

4.5 

0.7625 

5.7 

0.6,8 

4.0 

0.9681 

5.0 

0.8502 

4.55 

0.763 

5.8 

0.5973 

4.1 

0.954 

5.1 

0.85 

4.6 

0.7625 

5.9 

0.5744 

4.2 

0.9364 

5.6 

0.85 

4.65 

0.7611 

6.0 

0.5467 

4.25 

0.9261 

5.8 

0.8574 

4.7 

0.7585 



4.3 

0.9154 

5.9 

0.8646 




4.4 

0.8949 

6.0 

0.8735 





4.5 

0.8768 




x : Axial Position 

r : Radial Position 

R £ : Radius of Cowl Lip 


The first results employing the internal flow computational algo- 
rithm in which shock waves are discretely fitted are for the design con- 
ditions of M ot = 3.5, zero centerbody translation, and zero incidence 
(a - 0°). At the design point, the bow shock wave intersects the cowl 
lip exactly at zero incidence. Since the flow field is axi symmetric 
at zero incidence, it can be computed using a two-dimensional method. 
Comparisons of the results obtained from the present analysis with 
those obtained from a two-dimensional method of characteristics scheme 
(35) for the zero incidence design point conditions are shown in 
Figures 12 and 13. In these figures, the static pressure P normalized 
by the free-stream stagnation pressure Pj w is plotted versus the axial 
position x normalized by the cowl lip radius R e . Pressure distributions 
are shown for both the centerbody and the cowl. The results obtained 
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by the two-dimensional method of character)* sties algorithm are indi- 
cated by solid lines, and the results obtained by the present analysis 
are indicated by the dashed lines. Fifty radial stations were used in 
the two-dimensional method of characteristics solution. Figure 12 
illustrates the ease where a total of 11 radial stations (9 streamline 
points and an upstream and downstream shock wave point) were employed 
in the three-dimensional method of characteristics solution. Good over- 
all agreement is observed, A slight smearing of the pressure distri- 
bution downstream of the second intersection of the shock wave with the 
centerbody and a slight shifting of the shock wave-solid body inter- 
sections are present in the three-dimensional algorithm's results. 

The smearing of the pressure distribution is primarily a consequence 
of the coarse mesh size used in the three-dimensional scheme's solution. 
Figure 13 illustrates the solution obtained by the three-dimensional 
analysis when a total of 21 radial stations were used in the computa- 
tion. In this case, the agreement between the three-dimensional 
analysis and the two-dimensional analysis is excellent. The pressure 
distribution behind the second shock wave-centerbody intersection is 
predicted very well. The axial locations of the shock wave-solid 
boundary intersections also agree very well. For this computation, 
the maximum deviation in the computed mass flow rate at any solution 
plane from that at the cowl lip solution plane was approximately 0,77 
percent. 

Comparisons of the results of the present analysis with experi- 
mental data (36) for the Boeing Mach 3.5 inlet for a = 0° are shown in 
Figure 14. Generally speaking, good agreement is observed. The three- 

dimensional method of characteristics scheme predicts shock wave-solid 
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boundary intersections slightly downstream of the locations where the 
experimental data indicate the intersections to occur. Since the 
presence of a boundary layer would move the predicted intersection 
points forward, this result seems plausible. Note that the best agree- 
ment is obtained away from the regions where the boundary layer is 
removed {see Figure 11). 

At a given free-stream Mach number, the centerbody assembly must 
be translated forward of its design point position as the angle of 
incidence is increased to maintain supersonic flow through the geo- 
metric throat of the annulus. The forward translation of the center- 
body causes the cross-sectional area of the geometric throat to in- 
crease. Moreover, as the free-stream Mach number is reduced from the 
design point value, even further forward translation of the centerbody 
is required. An experimentally obtained centerbody translation 
schedule (37) is presented in Figure IB, where the nondimensional 
centerbody translation is denoted by Ax/R c - The effects of an increase 
in the angle of incidence and a reduction of the free-stream Mach num- 
ber are clearly illustrated in this figure. 

Results are presented below for two off-design conditions: 

(1) M^ - 2.5 with a centerbody translation of Ax/R e = 0.855, and 

(2) M - 3.3 with a centerbody translation of Ax/R„ = 0.356. For each 
of these off-design conditions, the internal flow field is computed 
for incidence angles of a = 0°, 3.0°, and 5.0°. For both off-design 
conditions, the results of the present analysis are compared with 
experimental data for an incidence angle of a = 3.0°. 

Results for the first off-design case of M^ - 2.5 and Ax/R c = 

0.855 are presented in Figures 16 to 19. Figure 16 illustrates the 
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computed centerbody and cowl pressure distributions for an incidence 
angle of a = 0°. Although the centerbody has been translated forward, 
the coordinate system origin is maintained at the forebody tip. Con- 
sequently, the internal flow computational regime begins at k/R = 
3.715. Generally speaking., the strength of the internal shock wave 
system for this case is somewhat reduced as compared to the design 
point case (see Figure 13). Figure 17 illustrates the computed 
pressure distributions and some experimental data for an incidence 
angle of a - 3,0°, Pressure distributions for the centeroody and the 
cowl on both the leeward and the windward meridians are shown. Com- 
pared to the a = 0° case, the strength of the internal shock wave 
system is increased on the leeward side but reduced on the windward 
side. Experimental data are presented for the centerbody pressure on 
the leeward meridian and for the cowl pressure on both the leeward and 
windward meridians. Generally speaking, good overall agreement 
between theory and experiment is obtained except in regions of high 
viscous interaction and boundary layer bleed. For all of the three- 
dimensional computations, 21 ci rcumferential stations and 11 radial 
stations (9 streamline points and an upstream and downstream shock 
wave point) were employed in the computed sector (half-plane). The 
maximum deviation of the mass flow rate at any solution plane from the 
mass flow rate at the cowl lip solution plane for the a = 3.0° case 
was 0.44 percent. The computed pressure distributions on the center- 
body and the cowl for both the leeward and windward meridians for the 
incidence angle of a = 5.0° are shown in Figure 18. The leeward 
meridian shock wave strength has been increased over the a = 3.0° case 
whereas the shock wave strength on the windward meridian has been 
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reduced. The maximum deviation in mass flow rate for the a = 5.0° case 
was 0.89 percent. Finally, to illustrate the effect of increasing 
angle of attack on the center body pressure distribution, the center- 
body results of Figures 16, 17, and 18 are superimposed in Figure 19. 

Results for the second off-design case of M m = 3.3 and Ax/R. = 
0.356 are presented in Figures 20 to 22. The computed pressure dis- 
tributions for the centerbody and the cowl for an incidence angle of 
a = 0° are presented in Figure 20. With the prescribed centerbody 
translation, the internal flow computational regime beings at x/R c = 
3.216. Figure 21 illustrates the computed centerbody and cowl static 
pressure distributions on both the leeward and windward meridians for 
an incidence angle of a = 3.0°. The strengthening of the leeward side 
shock wave and the weakening of the windward side shock wave are again 
noted. Experimental data for the leeward meridian of the centerbody 
and for both the leeward and windward meridians of the cowl are also 
shown in Figure 21. Fairly good overall agreement between theory and 
experiment is obtained until regions of high viscous interaction and 
boundary layer removal are reached. Again, 21 circumferential sta- 
tions and 11 radial stations were used in the computation. The maximum 
deviation of the mass flow rate at any solution plane compared to that 
on the cowl lip solution plane for the a - 3.0° case was 0.67 percent. 
Figure 22 illustrates the computed static pressure distributions for the 
centerbody and the cowl for an incidence angle of a = 5,0°. The ma:i- 
mum deviation in mass flow rate for this case was 0.84 percent. 

Finally, comparisons are made between the results of the present 
analysis and results obtained from the finite difference s frock- 

4 

capturing algorithm developed by Presley (37). At present, the 
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computer program developed by Presley is the only other analysis capable 
of predicting the internal flow field in supersonic mixed-compression 
inlets at angle of attack. That algorithm employs the second-order 
accurate finite difference operator devised by MacCormack (38). In 
that scheme, shock waves are automatically captured in the computa- 
tional mesh without requiring any special logic which discretely fits 
discontinuities. The presence of shock waves in the solution is evi- 
denced by steep gradients in the computed flow property fields. 

Figure 23 compares the centerbody and cowl pressure distributions 

obtained by the method of characteristics scheme to those calculated 

* 

by the shock-capturing technique for the case of = 3.3, a = 3.0°, 
and Ax/R = 0.356. For the most part, good agreement between the two 

V 

analyses is obtained. In the method of characteristics solution, 
however, the shock wave solid boundary intersections are more sharply 
defined. This result is to be expected, since in the shock-capturing 
technique shock waves are not discretely fitted but rather are smeared 
over a number of mesh points. The shock- capturing algorithm employed 
11 circumferential stations and 21 radial stations for the solution 
presented in Figure 23. 
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SECTION VII 


CONCLUSIONS 


The flow field in a supersonic mixed-compression aircarft inlet 
at nonzero angle of attack has been computed using the method of char- 
acteristics for steady three-dimensional flow in conjunction with a 
discrete shock wave fitting p :>cedure. The culmination of the present 
research effort is a production type computer program which has the 
capability to predict the flow field in a variety of axi symmetric 
mixed-compression aircraft inlets. A number of conclusions concerning 
the present analysis can be made: 

1. The external flow field about the forebody can be accurately 
calculated if a bow shock wave of reasonably strong strength 
exists. 

2. For axi symmetri c flows, the solution obtained by the present 
analysis agrees well with the solution obtained by the 
two-dimensional method of characteristics. 

3. Except in the regions of strong viscous interaction and 
boundary layer removal, the results of the present analysis 
agree well with experimental data. 

4. Good agreement is obtained between the present analysis and 
a finite difference shock-capturing technique for three- 
dimensional flow solutions. The present analysis, however, 
which discretely fits shock waves, provides a better resolu- 
tion of the shock wave structure. 
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5. Without the matching of the present analysis to a higher- 
order boundary layer analysis, including the influence of 
molecular transport in the computation has little or no 
effect on the solution. 

Although the inlets analyzed were axisymmetric inlets, the com- 
puter program can be readily modified to analyze geometries which have 
noncircular cross-sections. Moreover, the inclusion of finite rate 
chemiral reactions in the thermodynamic model is reasonably straight- 
forward. The analysis can be modified to compute the external flow 
about a stepped cone and to compute the internal flow when the bow 
shock wave has been ingested into the annulus. Perhaps the most worth- 
while endeavor, though, would be to mate the present analysis with a 
three-dimensional compressible turbulent boundary layer analysis. The 
boundary layer analysis should have well developed three-dimensional 
turbulence models, an accurate means of computing an oblique shock 
wave-boundary layer interaction in three-dimensions, and the capability 
to account for boundary layer removal. 



90 


i : 

^ Si .irjfsaratt. - *&>«&& ? 


liSiZX-t :£a£g£tait , .MXatwc:&* -V- -t4asi>4aj. jc:-ai;»«*. c&i'ij’ut -*gu&gtM. -■ 


APPENDIX A 
GOVERN IMG EQUATIONS 


1 . INTRODUCTION 

The major assumptions constituting the gas dynamic model are: 

1. continuum flow 

2. steady flow 

3. negligible body forces 

4. thermi.dynamic equilibrium (i.e., mechanical, thermal, and 
chemical equilibrium) 

5. no mass diffusion 

6. negligible radiative heat transfer and no internal heat 
generation other than viscous dissipation 

7. viscous and thermal diffusion effects of secondary impor- 
tance in determining the solution 

The governing equations for the assumed flow model consist of the con- 
tinuity equation, the component momentum equations, the energy equation, 
the thermal and caloric equations of state, and the appropriate 
representations for the molecular transport properties. These rela- 
tions are presented in this appendix 

2. DIFFERENTIAL EQUATIONS OF MOTION 

★ 

The general continuity equation (29) is 

y — - ' 

Repeated indices imply summation over the range of 1 to 3 unless other- 
wise noted. 
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Dp 

Dt 


+ p 



« 0 


(A.1) 


where t denotes time, p is the density, (i=l,2,3) denotes the three 
rectangular coordinates x, y, and z, respectively, and u.. (i=l,2,3) 
denotes the corresponding velocity components u, v, and w, respectively. 
The operator D( )/Dt in equation (A.l) is the material derivative 
given by 



(A.2) 


For steady three-dimensional flow, equation (A.l) may be written in 


expanded form as 


pu x + pv y + pw z + up x + vp y + wp z =0 (A. 3) 

where the subscripts x, y, and z denote partial different! on with re- 
spect to the corresponding direction. 

~ r he appropriate momentum equation is the Navier-Stokes equation 
(29), which written in component form is 


p 


Du. 


Dt 


1 - B, 





' 3U i . 

“ 

2 a 

r au : 

u 

ax. ax. 


" 3 3x. 



, J <j 

J 

1 

l 3J 


+ 




0 = 1 , 2 , 3 ) 


(A. 4) 


where denotes the ith component of the body force, Pis the pressure, 
y denotes the dynamic viscosity, and rj is the second coefficient of 
viscosity. 
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A major assumption of the present analysis is that the effects 
of viscous and thermal diffusion are of secondary importance in de- 
termining the solution as compared to the inertial effects. Consistent 
with this assumption of inertial dominance, the viscous and thermal 
diffusion terms in the governing differential equations will be treated 
as forcing or correction terms in the method of characteristic scheme 
to be presented. In the following, the viscous and thermal transport 
terms will be placed on the right-hand sides of the respective govern- 
ing equations. The convective terms will be placed on the left-hand 
sides, and will be considered as constituting the principal parts of 
these equations. Thus, writing equation (A. 4} with the assumptions of 
steady flow, negligible body forces, n ~ 0 [Stokes's hypothesis (30)]» 
and inertial dominance gives 


pu 


3U 1 + 3P 


j 3X 


3x. 


F i 


(1=1 ,2.3) 


(A. 5) 


where 

F. = — 
r i 3x . 

w 

Treating the viscosity as a variable, equations (A. 5) and {A. 6) 
can be written in expanded form for each of the three coordinate 
directions as 



3u. 3U. 

L 4 J. 


_ 2 a 

3U. 

1 1 




i 3x j N 


3 Sx. 



1 

Jj 


(1-1 .2,3) (A. 6) 


puu x + pvUy + pwu z + P x = (A. 7) 
puv x + pvVy + pwv z + P y = F y (A. 8) 
puw x + pvw y + pww z + P z = F z (A. 9) 
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where 


f x‘ v x 


n u .1 

3 x 3' ’y 


- -Kv.. + w 2 ) + p y (u y + v x ) + y z (u z + w x ) 


+ w {j “xx + “yy + U Z 2 + I< v ,y + w xzj] 

F y * “yl v y ' K + w z>] * >*x {v x + “y> + l 'z <v z + V 


(A. TO) 


+ V 


■k v + V + V + 
3 yy xx zz 


I<V + V>_ 


(A.H) 


F = u 
z ^z 


J w z - f^x + v y>] + “x (w x + u z> + “y (w y + V 


1 w zz + "w + w yy + T (u zx * v zy 


“I 

) 

— 


(A. 12) 


Finally, it remains to obtain an appropriate form of the energy 
equation. It is assumed in the present analysis that the working gas 
am be represented as a simple system in thermodynamic equilibrium. 
Under this assumption the thermodynamic relation (31) 


Tds = dh - 


dP 


(A. 13) 


is valid, where T denotes the absolute temperature, s is the entropy 
per unit mass, and h is the enthalpy per unit mass. For a simple 
system, specification of any two independent thermodynamic properties 
defines the thermodynamic state of the system (31). Thus, 


P = P(p,s) 


(A. 14) 


Employing the concept of the total derivative, and introducing the ma- 
terial derivative operator given by equation (A, 2), the following 
relation may be obtained from equation (A. 14). 


Ds 

Ot 


(A. 15) 


DP . fapl Dp . fSPl 
Dt - l^J s Dt * l3?J 

The sonic speed a is defined 


by 



Thus, equation (A. 15) may be written as 


(A. 16) 


DP 



(A, 17) 


where 



_3P) Ds 

3s Dt 

P 


(A. 18) 


The material derivative of entropy in equation (A. 18) may be 
expressed in terms of a thermal conduction function and a viscous dissi- 
pation function. Consider the energy equation in the following form 


( 29 ). 




POP 

Dt 


(A. 19) 


In equation (A. 19), e denotes the internal energy per unit mass, k is 
the thermal conductivity, and represents the viscous dissipation func- 
tion which for n = 0 is given by 


* * 


3u. 3U 

— 1 + 


2 > u k 




6 . . 


3x- 3 3x k ij 


(A. 20) 


where 6.. is the Kronecker delta. Using the definition of enthalpy 

* J 

(h = e + P/p) in equation (A. 13) yields 
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(A. 21) 


D 

Tds = de — £ dp 

P 

From equation (A. 21) the material derivative of internal energy may be 
written as 


Oe _ T Os , P Op 
Dt’ f 0t 2 Dt 


(A. 22) 


Introducing equation (A. 22) into equation (A. 19) yields 


Ds _ _3_ 
pT Dt * 


iT_ 

3X. 


+ $ 


Substituting equation (A. 23) into equation (A. 18) gives 


(A. 23) 




d 

3T 


3x* 

K 3x. 

+ $> 


(A. 24) 


where 


K = 


1_ 

pT 



(A. 25) 


By treating the thermal conductivity as a variable, and assuming 
steady three-dimensional flow, equations (A. 17) and (A. 24) may be 
written as 


UP X + vP y + WP 2 - a 2 {up x + v Py + WP 2 ) 55 F e (A. 26) 

where 

F e " ^{ K ^xx + T yy + T zz^ + K x T x + *y T y + *z T z 

+ m[*( 4 ♦ »y + + V, + V* * V,) + v x + W x 

t(u + v + w,)*U (A. 27) 

y y z z 3' x y z J J 


As in the component momentum equations, the viscous and thermal 
diffusion terms in the energy equation have been placed on the right- 
hand side and will be treated as forcing functions in the method of 
characteristics scheme to be presented. The left-hand side Is composed 
of the convective terms which are considered to constitute the princi- 
pal part of this equation. 

3. THERMODYNAMIC MODEL 

Before a solution to the system of governing partial differential 
equations can be obtained, the temperature T, sonic sp ad a, thermo- 
dynamic parameter £, viscosity y, and thermal conductivity *c must be 
expressed in terms of the dependent variables P and p. The representa- 
tions for T, a, and £ are discussed in this section. The relations 
for y and k are presented in the next section. 

The general functional forms of the temperature T, sonic speed a, 
and thermodynamic parameter £ may be expressed as 


T - T(P,p) 

(A. 28) 

a * a(P,p) 

{A. 29) 

e - c(p.p) 

{A. 30) 


For multicomponent systems, with either frozen or equilibrium chemical 
composition, the functional relationships for T, a, and £ are obtained 
from thermochemical calculations. In the case of a thermally and 
calorically perfect gas, the functional relationships for T, a, and £ 
are simple analytical expressions given by 
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T = P/pR 


(A. 31) 


a s (Tf/P) 1/Z " (A. 32) 

€ = y - 1 x. (A. 33) 

where y Is the specific heat ratio, and R is the gas constant, 

In the computer program developed in the present investigation, 
the temperature, sonic speed, and thermodynamic parameter £ are calcu- 
lated in a separate subroutine. The assumed thermodynamic model is 
that of a thermally and calorically perfect gas, thus, equations (A. 31) 
to (A. 33) are employed. Substitution of a replacement subroutine for 
the existing one allows other thermodynamic models to be specified. 


4. TRANSPORT PROPERTIES 

Representations are required for the viscosity, the thermal 
conductivity, and their spatial gradients. Both viscosity and thermal 
conductivity are functions of temperature and pressure. Hence, 

P = i«|T,P) (A. 34) 

k = k{T,P) (A. 35) 


Using equations (A. 34) and (A, 35), the spatial derivatives of viscosity 
and thermal conductivity may be written as 



(A. 36) 
(A. 37) 


98 


Hence, spatial derivatives of pressure and temperature are also re- 
quired. Spatial derivatives of pressure and density are employed In 
the basic integration scheme {even for the invlscid flow case). Thus, 
those derivatives are already available. Spatial derivatives of 
temperature can be expressed in terms of spatial derivatives of pressure 
and density by differentiating the thermal equation of state, equa- 
tion (A. 28). 

The pressure dependency Indicated in equations (A. 34) and (A. 35) 
is usually quite weak, and often both the viscosity and the thermal 
conductivity are assumed to be functions of temperature only. Thus, 


l> = v(T) 
K = <(T) 


(A, 38) 
(A. 39) 


The Sutherland formula (30) is a good representation for equation 
(A. 38). 


Vi = Vi, 


f T 1 

1.5 

T + S 
0 

T 


T + S 

oj 


h 4 


(A. 40) 


In equation (A. 40), vj q is the viscosity at the reference temperature 
T q , and S is a constant. Equation (A. 39) can be represented by the 
quadratic expression 


< - a-j + agT + a^T 


(A. 41 ) 


where the coefficients a^ (i=l,2,3) are obtained by curve fitting 
thermal conductivity data. 

In the computer program, the viscosity, the thermal conductivity, 
and their spatial derivatives are computed in a separate subroutine. 
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The assumed functional forms of viscosity and thermal conductivity are 
given by equations (A. 40) and (A. 41), respectively. The coefficients 
a^ (i=l,2 r 3) in equation {A. 41) are internally generated in the computer 
program by curve fitting user supplied data. Different formulations 
for the transport properties can be implemented into the computer pro- 
gram by supplying an appropriate replacement subroutine. 
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APPENDIX B 


DERIVATION OF THE EQUATIONS FOR THE CHARACTERISTIC 
SURFACES AND THE COMPATIBILITY RELATIONS 


1. INTRODUCTION 

Systems of hyperbolic partial differential equations in n inde- 
pendent variables have the property that there exist surfaces in 
n-space on which linear combinations of the original differential 
equations can be formed that contain derivatives only in the surfaces 
themselves. Differentiation in these surfaces is performed in (n-1)- 
space. The resulting differential operators are interior operators 
which are known as compatibility relations. The surfaces are called 
characteristic surfaces. A compatibility relation is valid only when 
it is applied on its corresponding characteristic surface. Furthermore 
data cannot be arbitrarily specified on a characteristic surface, but 
instead must satisfy the compatibility relation. 

The method of characteristics is based on replacing the original 
system of partial differential equations with an equivalent number of 
compatibility relations applied on the appropriate characteristic sur- 
faces. In flows with two independent variables, the method of char- 
acteristics has the advantage of reducing the solution of a system of 
partial differential equations to the solution of a system of ordinary 
differential equations. In three-dimensional flow, however, the 
resulting compatibility relations are still partial differential equa- 
tions in two independent directions. 


In this appendix, the equations for the characteristic surfaces 
and the corresponding compatibility relations are derived for steady 
three-dimensional flow. For a complete discussion of hyperbolic partial 
differential equations in three independent variables, reter to Courant 
and Hilbert (39). An excellent presentation of the method of character- 
istics for three-dimensional flow is given in Zucrow and Hoffman (5 ). 

2. EQUATIONS OF MOTION 

The partial differential equations of motion for steady three- 
dimensional flow consist of the three component momentum equations, 
the continuity equation, and the energy equation. Those equations are 
developed in Appendix A, and are repeated below for reference. 


puu x + pvu y + pwu z + P x = F x 

(B.l) 

puv x + pvv + pwv z + P y * F y 

(B.2) 

puw x + pvw y + pww z + P z = F z 

(B.3) 

pu x + pv y + pW z + up x + Vp y + wp z = 0 

(B.4) 

2 

uP x + vP y + wP z - a (up x + vp y + wp z ) - F g 

(B.5) 


In equations (B.l) to (B.5), u, v, and w denote the x, y, and z compon- 
ents of velocity, respectively, p is the density, P is the pressure, 
a is the sonic speed, and the subscripts x, y, and z denote partial 
differentiation in the corresponding direction. The nonhomogeneous 

terms F , F , F , and F are the forcing terms in the x, y, and z com- 
x y z 0 

ponent momentum equations and the energy equation, respectively. Writ- 
ten in this form, with the left-hand sides constituting the prinicpal 
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parts, equations (B. 1 ) to (B.5) may be classified as a system of quasi- 
linear nonhomogeneous partial differential equations of first order. 

The system is hyperbolic (i.e., has real characteristic surfaces) if the 
flow is supersonic. 


3. CHARACTERISTIC SURFACES 

The general compatibility relation, which is a linear combination 
of the governing partial differential equations, is formed by multiply- 
ing equations (B.l) to (B.S) by the arbitrary variables w. (i=l to 5), 
respectively, and summing. This yields 


<d 1 (puu x + pvUy + pwu z + P x ) + u) 2 (puv x + pWy + pwv z + P ) 

+ o> 3 {puw x + pvWy + pww z + P z ) + ta>4(pu x + pv y + pw z 

-f Up + Vp + Wp ) + W,-[uP + vP + wP 
M x *y M z' 5 L x y z 

2 

- a (up x + vpy + wp 2 )] = m^F x + u) 2 F y + co 3 F z + W gF e (B.6) 


Equation (B.6) may be written as 

+ pvw-jU + + pU0J 2 V x + P^ Vui 2 + ^4^ 

+ P ww 2 v z + P^ 3 w x + pvt^w + p( ww 3 + u 4^ w z 

+ (w 1 + uw 5 )P x + (u 2 + Va, 5)Py + Ug + Wtog)P 2 

2 2 2 
+ utu)^ - a w 5 )p x + v(ui^ - a Wgjp^ + w{w 4 - a Wg)p z 

= -,F X + ^ + “ 3 f z + "5 F e 


y 


(B.7) 
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By noting the coefficients of the partial derivatives in equation 
(B.7), the following vectors may be defined. 


M, - 

[p(uw-j 

+ 01 4 ), 

P Vto-| , 

f 

Q. 

(B.8) 

V 

[puu>2* 

p(vw 2 

+ w 4 ). 

pWujg] 

(B.9) 

' 

[pUD^, 

pVw^» 

p(ww 3 

+ u> 4 )] 

(B. 10) 

V 

[(w 1 + 

"“5 1 * 

( “2 + 

VWg), (u > 3 + WUJg )3 

(B.ll) 

V 

[uU 4 ■ 

• <1 <!).) 

, v(w 4 

2 2 
- a o>g) , w(to 4 - a cog)] 

(B.12) 


The directional derivative of a function f in some direction 
4 = ^x**y ,J ^ 1S given by 


df=r 31 + o 

dt x 3x y ay z 3z 


(B.13) 


By considering equations (B.8) to (8.13), equation (B.7) may be written 
as 


du , dv . dw , dP .do 
dWj dW 2 dW 3 dW 4 dW^ 


“l F x * 


tooF + 

2 y 


w 3 F z 


“5 F e 


(B.14) 


where du/dW-j is the directional derivative of u in the direction, 
etc. 

On a characteristic surface, equation (B.14) reduces to an interior 
operator, that is, differentiation takes place in the surface itself. 

For this to occur, the vectors W.. (i=l to 5) must all lie in the 
elemental plane which is tangent to the characteristic surface at the 
point in consideration. This means that the vectors WL (i=\ to 5) are 
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linearly dependent. Let the normal to the characteristic surface be 


denoted by N * ^ x ’^y T ^z^' Hence, on the characteristic surface 


N • W, = 0 (i=l to 5) 


(B. 1 5) 


Equation (B.15) yields five linear homogeneous equations which may be 
written in matrix form as follows 


pU 

0 

0 

pN 

X 

0 


U)-| 

0 

pU 

0 

P N y 

0 


*2 

0 

0 

pU 

PN Z 

0 


*3 

N 

N 

N 

0 

u 


U). 

X 

y 

z 




4 





2 



0 

0 

0 

u 

-a U 


_“s 


= 0 


(B. 16) 


where 


U = uN + vN + wN 
x y z 


(8.17) 


Since the system given by equation (6.16) is homogeneous, a nontrivial 
solution exits only if the coefficient matrix is singular, which means 
its determinant must be zero. Evaluating the determinant and equating 
it to zero yields 


(pU) 3 [U 2 - a 2 {N 2 + N 2 + N 2 )] = 0 


(B. 18) 


Equation ( B . 18) is the characteristic equation for the original system 
of equations, equations (B.l) to (B.5). The form of equation (B.18) 
is that of a repeated linear factor and a quadratic factor. 
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Equating the two factors in equation (B. 18) to zero yields the 
equations of two real nonintersecting cones formed by the envelope of the 
characteristic normals at a point. Setting the linear factor in equation 
(B.18) to zero gives (the case of p = 0 is immediately dismissed) 

uN + vN + wN = 0 (B.19) 

a jr Z 

Equation (B.19) represents a degenerate cone formed by the envelope of 
characteristic normals at a point, each normal being orthogonal to the 
local velocity vector. Hence, equation (B.19) represents a plane 
normal to a streamline. The characteristic surface is the reciprocal 
cone to this degenerate cone of normals, and, hence, is also degenerate, 
consisting of line segments tangent to the streamlines. Characteristic 
surfaces with normal components satisfying equation (B.19) are called 
stream surfaces. The envelope of all stream surfaces at a point is a 
single pencil of planes whose axis is a streamline. A streamline may 
be represented by the following equations 

dx/dt = u dy/dt = v dz/dt = w (B.20) 

where t is the time of travel of a fluid particle along the streamline. 

Equating the quadratic factor in equation (B.18) to zero gives 

(uN x + vNy + wN z ) 2 - a 2 (N 2 + N 2 + N 2 ) = 0 (B.21) 

Equation (B.21) represents the quadric surface of a right circular cone 
formed by the envelope of characteristic normals at a point. In gas 
dynamics this cone is usually referred to as the cone of normals, and 
is a real cone if q > a, where q is the velocity magnitude. Equation 
(B.21) may be written as 
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a 


(B.22) 


un + vn + wn = 
x y z 

where n = (n ,n ,n ) is the unit normal to the characteristic surface. 
x x y 2 

Equation (B.22) was obtained by arbitrarily selecting the positive root, 
and the results which follow are consistent with that selection. Char- 
acteristic surfaces whose normal components satisfy equation (B.21), 
or equation (B.22), are called wave surfaces. 

Equation (B.21) is the equation for the cone of normals, which is 

* 

a quadric surface. In general, a quadric surface may be expressed as 
A. .dx.dx. = 0 (B.23) 

l j l j 

where x. ( i =1 ,2.3) denotes the three cartesian coordinates x, y, and 
2 , respectively, and A is a nine element coefficient matrix of order 
two. A normal vector is a directed line segment, so 

N. = o dx. { i=l ,2,3) (B.24) 

where N,. is the itfi component of the normal vector, and o is a constant 
proportional to the length of the normal. By considering equations 
(B.23) and (B.24), equation (B.21) may be written as 

(u.u. - a^. Jdx.dx. - 0 (B.25) 

i J i j i J 

where u.. (i = 1 ,2*3) denotes the three velocity components u, v, and w, 
respectively, and 6^ is the Kronecker delta. 


^Repeated indices imply summation over the range of 1 to 3 unless other- 
wise noted. 


107 


The characteristic cone, which is the envelope of all wave 
surfaces at a point, is the reciprocal cone to the cone of normals 
given by equation (B.21), or equation (B.25), The geometrical rela- 
tionship between these surfaces is shown in Figure B.l. If the general 
form of the equation of the cone of normals is given by equation 
(B.23), then the reciprocal cone is given by { 9) 

ATldx.dx. =0 (B, 26) 

U i J 

where is the inverse of the nine element symmetric matrix A in 
equation (B.23). Using equation (B.25) to determine A from wijich A~^ 
may be determined, equation (B.26) for the characteristic cone 
wri tten as 

[u.Uj - (q 2 - a Z )6 ij ]dx.dx j = 0 (B. 27) 

Equation (B.27) represents a real cone if q > a. Writing equation 
(B.27) in expanded form yields 

[u 2 - <q 2 - a 2 )]dx 2 + [v 2 - (q 2 - a 2 )]dy 2 + [w 2 - (q 2 - a 2 )]dz 

% 

+ 2uv(dx)(dy) + 2uw(dx)(dz) + 2vw{dy)(dz) = 0 (B.28) 

The characteristic cone given by equation (B.28) is known as the Mach 
cone and represents the envelope of all wave surfaces at a point. The 
line of tangency between a particular wave surface and the Mach cone 
is known as a bi characteristic. Integration of equation (B.28) gives 
the curved cone known as the Mach conoid. 

In summary, for steady three-dimensional flow there are two 
families of characteristic surfaces: stream surfaces and wave surfaces 



FIGURE B.l. RELATION BETWEEN THE CONE OF NORMALS AND 

THE CHARACTERISTIC CONE (MACH CONE) 


{see Figure B.2). The normal to a stream surface must satisfy equation 
(B.19), and, hence, the stream surface contains the local velocity 
vector. The envelope of all stream surfaces at a point is the streamline 
through the point. The normal to a wave surface must satisfy equation 
(B.21). The envelope of all wave surfaces at a point is the Mach cone. 
The line of contact between a particular wave surface and the Mach cone 
is called a bi characteristic. At any point there are an infinite 
number of stream surfaces and wave surfaces. 

4. SOLUTION FOR THE w. 

On a characteristic surface, equation (B. 1 4) reduces to an interior 
operator, that is, it becomes a compatibility relation. To obtain the 
exact form of the compatibility relation, the ^ (i=l to 5) must be 
determi ned. 

For a stream surface, equation (B.19), repeated below, is valid. 

uN + vN + wN_ = U = 0 (B.19) 

x y z 

Substitution of equation (B.19) into the homogeneous system given by 
equation (B.16) yields 
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0 

0 
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The coefficient matrix in equation (B.29) is rank two (rank is the 
number of nonzero rows in the row echelon form of a matrix). The 
number of independent nontrivial solutions for the u. is equal to the 
order of the coefficient matrix minus its rank, and hence, in this 
case, is three. From equation (B.29), = 0 for all solutions, w,. 

is arbitrary, while satisfy the following equation. 


“l N x 

+ 

V 

“3 N z ’ 0 

(B.30) 

three 

possible 

solutions is 


“1 = 

“2 = 

w 3 1 

= 0J 4 = 0, w 5 = 1 

(B.31) 

“1 = 

u. 

w 2 5 

s v, = w, u> 4 = u)g = 0 

(B.32) 

= 

V 

u> 2 

= Sy, ^3 = S 2 * - 0 

(B.33) 


The vector S = (S »S ,S ) in equation (B.33) lies in the stream surface 
x y z 

and is independent of the velocity vector. 

On a wave surface, equation (B.21) is valid. That equation may be 
written as 


U = a j N j (B.34) 

where |ff| Is the magnitude of the normal to the wave surface. Substi- 
tuting equation (B.34) into equation (B.T6) yields 


112 


pa | Ni| 0 0 pN 


0 pajNj 0 pN 0 a) 2 


0 0 pajNj pN z 0 


(B.35) 


N x N y N z G a l 


j^O 0 0 a | N | -a | Ny K 


The coefficient matrix in equation (B.35) is rank four, and, hence, one 
independent nontrivial solution exists for the w,.. The solutions for 
w.|, ujg* b>£* and Wg may be expressed in terms of od^. Arbitrarily 
selecting = -1 yields 

- n^/a, Wg - n y /a* = — “1* 


w 5 = -1/a 


(B.36) 


where n - (n x ,n y ,n z ) * s thG un ^ norma ^ tG wave surface. 

5. COMPATIBILITY RELATIONS 

The compatibility relations are obtained by substituting the 
solutions for the ^ into equation (B.6). The compatibility relations 
valid along the stream surfaces are obtained by substituting equations 
(B.31) to (B.33) into equation (B.6). The results are 


UP X + vp y + WP 2 - a (up x + v Py + wp z ) = F e 
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pu(uu + vu + wu ) + pv(uv + vv + wv ) + pw(uw u + vw 
A y * xyz x y 

+ ww ) + Up + vP + wP - uF + vF + wF (B.38) 

C A Jr 4 x y z 

pS x (uu x + vu y + wu z ) + ps y (uv x + vv y + wv 2 ) + ps z (uw x + vw y 

+ ww ) + $ P +SP +SP = S F +SF +SF (B.39) 
z xx y y z z xx yy zz ' ' 


Note that equation {B- 37) is the same as equation (B.5), which shows 
that the energy equation is characteristic to begin with. 

Equations (B.37) and (B.38) may be written in a form that repre- 
sents differentiation in the streamline direction only. From equation 

(B.13)» noting that for a streamline l - u, l = v, and z - w, 

x y z 

the directional derivative along a streamline is given by 

-SlLJ, = u JL + v i- + w — ^ (B 40) 

dt 3x 3y w 3z 

where t is the time of travel of a fluid particle along the streamline. 
Using equation (B.40), equations (B.37) and (B.38) may be rewritten as 


dP - a 2 &. 

dt dt 


(B.41) 


pu at + pv av + DW a? + ft = uF x + vF y + wF z (B - 42) 

The compatibility equation that is valid along wave surfaces is 
obtained by substituting equation (B.36) into equation (B.6). The 
result is 
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(B.43) 


pan (uu + vu + wu ) + pan (uv + vv + wv ) 

AAyz y x y z 

+ pan z (uw x + vw y + ww z ) + (an x - u)P x + (an y - v)P y 
+ {an^ - w)P z - pa Z (u x + v y + w 2 ) = * 

where 

A * a <"x F x + Vy + "z F z> - F e < B ' 44 > 

Equation (B.43) may be written in a form that contains differen- 
tiation in the bi characteristic direction. A bi characteristic is a 
ray or generator of the Mach cone. The Mach cone is the reciprocal cone 
to the cone of normals (see Figure B.l). As a consequence, a bichar- 
acteristic is orthogonal to the surface of the cone of normals. The 
equation for the cone of normals is given by equation (B.21). Substi- 
tution of equation (B.24) into equation (B.21) yields the equation for 
the surface of the cone of normals in standard form [f(x,y,z) * constant]. 
Differentiation of this expression to obtain the gradient yields the 
direction of the bi characteristic. This gives l ~ (u - an ), 

X X 

l - (v - an^), and t z = (w - an z ) in equation (B.13), so that differ- 
entiation in the bi characteristic direction is given by 

3p = < u * an x ) fp + (» - *n y ) |p + (« - an 2 ) fP (B.45) 

In equation (B.45), t is the time of travel of a fluid particle along 
the streamline that is the axis of the Mach cone. The relationship be- 

••v 

tween the vectors l, V, and n is shown in Figure B.3. 
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The term 


2 2 2 2 

±{pa [n u + n v + n w + (u + v }n n + (u + w )n n 
L x x y y z z ' y x xy ' z x' x z 

+ (v 2 + w y )n y n z ]} 

may be added to and subtracted from equation (B.43), and then by employ 
ing equation (B.45) the following form of the wave surface compati- 
bility relation may be obtained. 


du . dv , 
pan x dt + pan y dt + 

+ {n y “ 1}v y + (n z " 

+ n y n 2 (v z + w y )] 


pan zl|-f = X - pa2[(n x - 1)u x 


l)w + n n (u + v ) + n n (u + w ) 
' z x y y x 7 x z z x / 


(B.46) 


The terms in brackets in equation (B.46) are known as cross derivatives 
and represent differentiation in the wave surface in a direction normal 
to the bi characteristic direction. 

Equations (B.29) and (B.35) determine the number of independent 
differential compatibility relations valid along a particular stream 
surface and a particular wave surface, respectively. At any point 
there exist an infinite number of stream surfaces and wave surfaces. 
However, the number of independent compatibility relations cannot exceed 
the number of independent equations of motion. Hence, it is necessary 
to determine which of the possible combinations of compatibility rela- 
tions are independent. Rusanov (32), using a proof in the space of 
characteristic normals, has shown that for steady three-dimensional 
isentropie flow two of the stream surface compatibility relations and 
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the single wave surface compatibility relation applied along three 
different wave surfaces form an independent set of characteristic 
equations. Rusanov’s results may be extended to the present problem 
since the principal parts of equations (B.l) to (B.5) are the same as 
those for isentropic flow. Thus, for the present problem, an inde- 
pendent set of compatibility equations consists of equations (8.41 ) 
and (B.42) applied along a streamline, and equation (B . 43) [or equa- 
tion (8.46)3 applied along three different wave surfaces. 

6. BUTLER’S PARAMETERIZATION OF TNE CHARACTERISTIC EQUATIONS 

The numerical algorithm that is employed in the present investiga- 
tion is based on a second-order scheme devised by O.S. Butler (24). 

This scheme has been used by Ransom, Hoffman, and Thompson ( 9) to 
compute isentropic steady three-dimensional nozzle flows, and by Cline 
and Hoffman (25) to compute chemically- reacting steady three-dimensional 
nozzle flows. 

In this section, Butler's parameterization of the characteristic 
equations is presented. The discussion below Is limited to the partic- 
ular application of Butler's method to the present problem. An excel- 
lent review of Butler's general method is given in Ransom, Hoffman, and 
Thompson (9). 

For Butler's scheme to be applicable, the characteristic determin- 
ant must be composed of a quadratic factor and a repeated linear factor. 
The determinant of the coefficient matrix in equation (B.16) is the 
characteristic determinant for the present problem, and by examination 
of equation (B . 18) it is seen that it is composed of the required 
factors. The quadratic factor corresponds to the wave surfaces. The 
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envelope of all wave surfaces at a point is the Mach cone. The line of 
tangency between a particular wave surface and the Mach cone is a bi- 
characteristic. The linear factor corresponds to the stream surfaces. 
The axis of the envelope of all stream surfaces at a point is a 
streamline. Butler's method assumes that for the linear factor, 
differentiation can be expressed soley along the axis of the envelope 


of the corresponding characteristic surfaces. Examination of equations 


{ B . 41 ) and (B.4?) demonstrates that this condition is applicable. 

As discussed in the first section of this appendix, if the system 
of governing partial differential equations has differentiation oc- 
curring in n-space, then differentiation in the characteristic surfaces 
occurs in {n-l)-space (i.e., differentiation is performed in a mani- 
fold of one lower dimension). As a result, for three-dimensional flow 
(n=3), the general form of a compatibility relation valid along a 
characteristic surface may be written as 

E ,0 u /3xi) + F Ou /Sx') = D (B. 47) 

v v i v v 2 

where the repeated index v implies summation over the range of 1 to 5, 

( i = 1 , 2 ) denotes two independent directions in the characteristic 
surface, u y (v=l to 5) denotes the dependent variables, and E , 

F y {v=l to 5), and D are general functions of x.! and u y . For stream 
surfaces, differentiation may be expressed solely in the streamline 
direction [see equations (B.41) and { B- 42 ) ] . Consequently, in the 
following, the discussion will be limited to the wave surfaces. 

For steady three-dimensional flow, Butler introduced the following 
parametric representation for a bicharacteristic. 


< 


f 


| 

f 

\ 

* 


} 


t 
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{ B , 48 ) 


dx^ = (u i + coucos© + e$.sin©)dt (i=l,2,3) 

In equation (B.48), (i=l,2,3) denotes the three cartesian coordinates 

x, y, and z, respectively, u. (i=l,2,3) denotes the corresponding 
velocity components u, v, and w, respectively, 0 is a parametric angle 
denoting a particular element of the Mach cone and has the range 
0 < 9 < 2 tt, t is the time of travel of a fluid particle along the 
streamline that is the axis of the Mach cone, and c is defined by 

c 2 = a 2 q 2 /(q 2 - a 2 } (B.49) 

where q is the velocity magnitude and a is the sonic speed. The vec- 
tors eij and $. are parametric unit vectors with ou, ^ , and 
u-/q (1=1 .2,3) forming an orthonormal set. A geometrical representation 
of this parameterization is given in Figure B.4, 

The direction specified by equation {B.48) lies in the wave surface 
and is in the bi characteristic direction. A direction in the wave 
surface and orthogonal to the bi characteristic direction may be written 
in parametric form as 

m = cB^eos© - cousin© (i-1,2,3) (B.50) 

Verification of the orthogonality of the directions given by equations 
(B.48) and (B.50) may be accomplished by forming the dot product 
(m.dx^) and using the orthonormality relations 
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FIGURE B.4. BICHARACTERISTIC PARAMETERIZATION 
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(B.51) 


Vi * Vi 


a .0 • 
1 P 1 


0 


Vi = Vi * Vi /q • 1 


By considering equation (B.47) and selecting x^ and x£ as the 
directions given by equations (B.48) and (B.5Q), respectively, the 
following form of the wave surface compatibility relation is obtained. 

A (u- + ca.cosQ + c0 .sine) ~ 

V 1 I | dX.j 

9u 

8 + C (cB.cos© - ca.sin©} r— - (B.52) 

VI i dXj 

In equation (B.52), A v , B, and C y are functions of ©* u y , and x^. 
Employing equation (B.13), and noting from equation (B.48) that along a 
bi character! stic 


+ ccx.cos© + cB^sinS ( i =1 ,2.3) (B , 53) 

equation (B.52) may be written as 

du au 

A v " B + c v ( c ^i cose " ca.sin® ) (B. 54) 

■ i 

where the operator d(. )/d£ represents the directional derivative along 
the bi characteristic. The general forms of the coefficients A , B, 
and G v are given by Butler as 

A v = A lv + A 2 v cos0 + A 3v sine (8-55) 

B = + Bpcose + B 3 sin9 (B.56) 
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G v = G lv + G 2o cos0 + G 3v S * n0 


(B.57) 


where the A^, B^, and (k=l,2,3 and v=l to 5) are independent of 0 
In addition to the parametric wave surface compatibility relation, 
given by equation (B.54), Butler also developed a noncharacteristic 
relation which is applied along a streamline. This noneharacteristic 
relation is used in the numerical scheme in conjunction with the wave 
surface compatibility relation applied along four different bi char- 
acteristics, and permits the formulation of three independent linear 


combinations of these five equations which do not contain cross deriva 
tives at the solution point. The cross derivative terms [see equation 
(B.46)] represent differentiation in the wave surface but in a direc- 
tion orthogonal to the bicharacteristic direction [i.e., differentia- 
tion in the direction given by equation (B.50)]. Butler presents the 


noncharacteristic relation in the form 


du 3 u 

A TvdF = 8 1 + * C 2v ce i ' C 3v ca i } 3x7 


(B. 58) 


where the operator d( )/dA represents the directional derivative along 
the streamline. The coefficients A-j , and (v*l to 5) in 

equation (B.58) are obtained by inspecting the form of equation (B.54) 
and then using equations (B.55), (B.56), and {B.57}. 

For the present problem, the actual form of the parametric wave 
surface compatibility relation, equation (8.54), may be obtained by 
substituting the appropriate parametric form of the wave surface unit 
normal into the compatibility relation, equation (B.43). The normal 
to the wave surface is also the normal to the Mach cone at a point 
common to both surfaces. The quadric surface of the Mach cone is 
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j i"i 


represented by equation (B.27), repeated below. 
[u^Uj - (q 2 - a 2 )<S.,.]dx.dx,. * 0 


(B.27) 


Substituting the parametric form for dx^, given by equation (B. 48) . 
into equation (B.27) yields 

[u.u. - (q 2 - a 2 )5. .](u. + ca.cose + c6.sin0)dx. - 0 (B.59) 

lj * J J J J * 

The ith component of the normal to this surface is 
N i = [u.u^. - (q 2 - a 2 )6 ij .](u j + choose 


+ cp.sine) (i=l .2,3) 

J 


(B.60) 


Employing the orthonormali ty conditions given by equation (B.51), equa- 
tion (B.60) may be written as 

Nj ^ a Z [u i - (q 2 /c) (a,.eos0 + B^sine)] ( i - 1 * 2 , 3 ) (B, 61} 

Dividing equation (B.61) by the magnitude of the normal JjNj = 

(N^N,j and using equation (B.51), the parametric form of the wave 
surface unit normal is obtained. 

n i = (a/c) (cu^/q 2 - a. cose - B^inQ) {i =1 *2,3) (B.62) 

Substituting equation (B.62) and the orthonormality relation 


Vj + <¥j + u i u j /q r 6 ij 


(B.63) 


into the wave surface compatibility relation, equation (B. 43) » gives 
the following parametric form of that equation 
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(B.64) 


dP 


du. 


+ pentose + e.sine)^ 1 * * - pc' (a.stne 


3u. 


- e. cose) (a. sine - 3,cosek— - 

* J J oXj 


where 


$ = - (c 2 /a Z )X 


(B.65) 


The operator d( )/dt in equation (B.64) denotes differentiation in the 
bi characteristic direction. 

It should be noted that the directional derivatives in equations 
(B.46) and (B.64) are not identical. - The directional derivative in 
equation (B.46) is based on equation (B.4S). Substitution of the 
parametric unit normal, given by equation (B.62), into equation (B.45) 
yields 


* (a 2 /c 2 )(u i - + ca.cosB + cB^sin6)|^-^ (B.66) 

The directional derivative in equation (B.64) is given by 

= (u.j + co^cos© + cosine )|£-^ (B.67) 

2 2 

Hence, the two expressions differ by the factor (a /c ). 

Finally, it remains to determine the actual form of the nonehar- 
acteristic relation, equation (B.58). Denote u^ (v=l to 5) and 
(i=l,2,3) in equations (B.54) and (B.58) by 

U-j = u, u 2 * v, u 3 = w, u 4 * P, Ug = P 

x 1 = x, x 2 = y, x 3 s z (B.68) 
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By Inspection of equation (B.64), and use of equations (B.68), (B.55), 
(B.56), and (B.57), the noncharacteristic relation is seen to be 


HP 9 on,* 

at = 0 - pc (a i a j + Mjfe 

J 

where 

a - (c 2 /a 2 )F e - (c 2 /q 2 )(uF x + vF y + wF z ) 


{B. 69) 


(B.70) 


The operator d( )/dt In equation (B. 69) denotes the directional deriva- 
tive along a streamline. 

In summary, Butler has developed a bi characteristic parameteriza- 
tion given by equation (B.48). The corresponding parametric form of 
the wave surface compatibility relation is given by equation (B.64). 
Butler also developed a noncharacteristic relation, given by equation 
(B.69), which is applied along a streamline. These relations, along 
with the stream surface compatibility relations, equations (B.41) 
and (B.42), constitute the system of compatibility relations. The use 
of this system of equations in the various unit processes is presented 
in Appendix E. 



APPENDIX C 


INTERPOLATION 


1 . INTRODUCTION 

In the course of computing the flow field, a number of situations 
arise which require interpolation. To this end, univariate, bivariate, 
and tri variate interpolation polynomials are employed in the numerical 
algorithm. These interpolation schemes are presented in this appendix. 

2. UNIVARIATE INTERPOLATION 

Univariate interpolation is required in geometry description, 
calculation of the transport forcing terms, and in determination of the 
properties along a space curve formed by the locus of shock wave solu- 
tion points. Applications to geometry description and transport term 
computation are discussed in Appendices D and G, respectively. The 
application to the determination of properties along a shock wave is 
discussed here. 

When a shock wave intersects either a solid boundary or a solution 
plane (a plane of constant x), a space curve is defined as illustrated 
in Figure C.l. Interpolated values of position, shock wave angle, and 
flow properties are required along this curve. For this purpose, the 
quadratic polynomial 

f ( 0 ) = a 1 + a 2 @ + a 3 © 2 (C.l) 
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FIGURE C. I . POINT STENCILS FOR UNIVARIATE INTERPOLATION 


is employed, where f(0) denotes a general function expressed in terms 
of the polar angle 0 given by 

e = tan" 1 (z/y) (G.2) 

where y and z are the coordinates of a point on the space curve. The 
coefficients a i {1=1, 2, 3) in equation (C.l) are determined by fitting 
this expression to three data points on the space curve, and, as a 
consequence, a system of three simultaneous linear equations must be 
solved for the coefficients a^ of each function representation. The 
solution to this system of equations is obtained using a Gaussian 
elimination method with complete pivoting (40). 

Figure C.l illustrates typical data point stencils used for de- 
termining coefficients in equation (C.l). The fit point array con- 
sists of a base point, which is the point closest to the position of 
the interpolated point, and the immediate neighbors of the base point. 

3. BIVARIATE INTERPOLATION 

Bivariate interpolation is required for property determination in 
a given solution plane (a plane of constant x). Two types of bivariate 
interpolation polynomials are employed in the numerical algorithm. 

They are a linear bivariate polynomial whose three coefficients are 
determined by fitting this expression to three data points, and a 
quadratic bivariate polynomial whose six coefficients are determined 
by a least squares fitting of nine data points. 

The linear bivariate polynomial is used in the single appli- 
cation when a streamline-shock wave intersection point is sufficiently 
close to the current solution plane so that an interior point unit 
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process on the downstream side of the shock wave is not performed. 

In that case the projection of the streamline onto the solution plane 
and subsequent property interpolation in this plane is performed. The 
bivariate interpolation polynomial used in this case is 

f{y»z) = + a^ + a 3 z (C.3) 

where f(y,z) denotes a general function of the coordinates y and z. 

The coefficients (i=l,2»3) in equation (C.3) are determined by 
fitting this expression to three data points. This yields a system 
of three simultaneous linear equations for the coefficients a,, of each 
function representation. This system of equations is solved using a 
Gaussian elimination method with complete pivoting [as was done for 
equation (C.l)]. 

A typical data point stencil used for determining the coefficients 
in equation (C.3) is illustrated in Figure C.2. Two shock wave solu- 
tion points and a field point constitute the fit point array. 

In all other situations which require bivariate interpolation, the 
quadratic polynomial 

2 2 

f(y,z) = a.j + a^y + a 3 z + a^yz + a g y + a g z (C.4) 

is employed, where f(y,z) is a general function of the coordinates y 
and z. The coefficients a.. (i=l to 6) in equation (C.4) are determined 
by a least squares fit of nine points. Using the standard theory of 
least squares (40), the system of normal equations which determines 
the coefficients in equation (C.4) is 
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FIGURE C . 2 . POINT STENCIL FOR LINEAR BIVARIATE 

INTERPOLATION 
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In equations (C.5) to (C.10), the £ sign implies summation over the 
range of 1 to 9, while the subscript i denotes the i th data point 
(i=l to 9). This system of simultaneous linear equations has a sym- 
metric coefficient matrix and is solved using a Gaussian elimination 
method with pivoting in the main diagonal. 
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Figure C.3 illustrates typical data point stencils used in de- 
termining the coefficients in equation (C.4). Basically, there are two 
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types of stencils: interior point and boundary point. Since the 
shock wave mathematically represents a discontinuity, the boundary 
point stencil must be employed when the interpolation base point (the 
data point closest to the interpolated point) is on the shock wave. 

The fit point array consists of the base point and its eight immediate 
neighbors. Special logic in the computer program is used to insure 
that no stencil bridges the shock wave. 

4. TRI VARIATE INTERPOLATION 

Trivariate interpolation is required for property determination 
on the surface of a solid boundary (a stream surface) and for property 
determination on the upstream and downstream sides of the shock wave. 
Two types of tri variate interpolation polynomials are employed in the 
numerical algorithm. They are a linear tri variate polynomial whose 
four coefficients are determined by fitting this expression to four 
data points, and a quadratic tri variate polynomial whose eight coef- 
ficients are determined by a least squares fitting of fourteen data 
points. 

The linear tri variate polynomial is used in the single applica- 
tion for property determination on the upstream side of the shock wave 
surface. This polynomial has the form 

f(x,y,z) = a 1 + agX + a^y + a^z (C.ll) 

where f{x,y,z) is a general function of the coordinates x, y, and z. 
The coefficients a^ (i=l,2,3,4) in equation (C.ll) are determined by 
fitting this expression to four data points. Hence, a system of f our 
simultaneous linear equations must be solved for the coefficient a i 
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of each function representation. This system of equations is solved 
using a Gaussian elimination method with complete pivoting [as was 
done for equations (C.l) and (C.3)]. 

A typical data point stencil used for determining the coefficients 
in equation (C.ll) is illustrated in Figure C.4. Three data points 
are located on one space curve and one data point is located on the 


other space curve. 

In all other situations which require trivariate interpolation. 


the quadratic polynomial 

f(x,y,z) = a ] + a^y + a g z + a 4 yz + a & y l 2 + a g z 2 


+ a^xy + a Q xz 


(C. 12) 


is employed, where f(x,y,z) is a general function dependent on the 
coordinates x, y, and 2 . The coefficients a. (i=l to 8) in equation 

nr 

(C.12) are determined by a least squares fit of fourteen data points. 
From the theory of least squares, the system of normal equations 
which determines the coefficients in equation (C.12) is 

14a, + l V 2 + 2 z i a 3 + l Vi a 4 + J y 1 a 5 + E z i a 6 


+ y x.y.a, + y x.z.a Q =7 f. 
L ri 1 L i i 8 L 1 


(C. 1 3) 


l y,a, + l yfa 2 + l y,z,a 3 + l y^z.a 4 ♦ J y?a 5 t £ y^a,. 


9 

+ y x.y.a. J r y x.y.z.a- = y y.f. 
L Y'l 7 L i t i 8 L i 


(C.14) 
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INTERPOLATION 


. 3&S9S3£^&'ia^^ ■ <#^i^‘-^ i &^-il , .jciiSrS« l 




'. .w_Wd ■ . -*-**&*$ „y .“r?-i rr-^.yir • -HC¥»J xX&twn ■ 


l z.a, + I y i z i a 2 + l 4> 3 + l y,.z-a 4 + I y,Vs + Z z i a e 

* l Wi a 7 + Z x i z f a 8 ’ l Z i f i (C ' ,S) 


Z ^i z 1 a T + l yWz + 1 y i Z i a 3 + l y i z i a 4 + l y i z 1 a S 

+ Z y f z ? a 6 + Z Vi z i a 7 + Z x t y i z i a 8 = Z y i z . f i <c - 16) 


Z yfa, + z yj a 2 + Z y^ 1 a 3 + Z y? z i a 4 + Z y? a s + Z y? z v> 6 
+ Z x,.y?a 7 - z Vi z t a 8 ■ Z y?f, (c.i?) 


Z z j a , + Z yi z i a 2 + Z z ? a 3 + Z y, z ? a 4 + Z y i z i a 5 + Z z ?6 


+ Z Wl a 7 * Z x i z i 8 = Z z ,' f i 


(C.18) 


z x ,y,«, + z Xjy? a 2 + z Vi z i a 3 + Z x i y i z i a 4 + Z Ws 


+ Z x,-yi z i a 6 + Z x ?yf a 7 + Z x^-Vs = Z x/if, (C.19) 


Z V,a, + 

+ Z Xj z j a 
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(C . 20) 


In equations (C. 1 3} to (C.20), the l sign implies summation over the 
range of 1 to 14, while the subscript i denotes the ith data point 
(i=l to 14). This system of simultaneous linear equations has a sym- 
metric coefficient matrix and is solved using a Gaussian elimination 
method with pivoting in the main diagonal [as was done for equation 
(C . 4 ) ] . 
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Figure G.5 illustrates typical data point stencils used in de- 
termining the coefficients in equation (C.12). The fit point array 
consists of seven data points along each of the appropriate space 

curves on either the shock wave or the solid boundary. 

It should be noted that a ten term quadratic trivariate poly- 
nomial, with coefficients determined by a least suqares fit of fourteen 
data points, was tried in place of equation (C.12). Use of this 
polynomial in flows with axial symmetry, however, did not produce 
results which were as symmetrical (especially for transverse velocity 
components) as those produced by equation (C.12). This result could 
possibly be due to the effects of ill-conditioning as is discussed in 
Hanming (40). Furthermore, scaling of the dependent variables did 
not appear to produce any improvement. 
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FIGURE C . 5 . POINT STENCILS FOR QUADRATIC TRIVARIATE 

INTERPOLATION 
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APPENDIX D 

SURFACE REPRESENTATIONS, AND STREAMLINE- AND 
BICHARACTERISTIC-SURFACE INTERSECTIONS 


1. INTRODUCTION 

The procedures employed for representing the solid boundary and 
shock wave surfaces are presented in this appendix. The technique used 
for determining the intersection point of either a streamline with the 
shock wave, or a bi characteristic with either the shock wave or the 
solid boundary, is also discussed. 

2. SOLID BOUNDARY SURFACES 

The centerbody and cowl surfaces are specified in the computer 
program by a separate geometry module that has the capability to de- 
scribe a variety of axi syimnetri c contours. More arbitrary geometries, 
such as those having elliptical or superelliptical cross sections, may 
be considered by supplying an appropriate replacement module. In 
general, to specify a surface completely, its functional form 
[f(x,y,z) = constant] and its gradient at any point [Vf(x,y,z)] must 
be available. 

The existing geometry module, which describes axi symmetric con- 
tours, divides the axial (x) domain into a number of intervals. In any 
interval, the body radius may be specified by either tabular input, 
cr by supplying the coefficients in a cubic polynomial written as a 


function of x. For the tabular input ease, linear interpolation is 
performed to obtain the radius r(x) between the points (x^,r^) and 
» r .| + ] ) where (x^ _< x < x. + ^). Alternatively, employing the 
cubic polynomial 

r(x) = a. + b.(x - x^ + c.(x - x^ 2 + d^x - x.) 3 

(Xj<x<x. +] ) (D.l) 

requires that the coefficients a.., b., c^ , and d^ be supplied for the 
ith interval (these coefficients must be externally generated). Since 
equation (D.l) is a cubic, slope and curvature can be matched at the 
junction point between two adjacent intervals (i.e., spline fits can 
be employed). 

3. SHOCK WAVE SURFACE 

Some of the unit processes, which are described in Appendix E, 
require an analytical representation for the shock wave surface. 

During the course of the program development, a number of different 
representations were devised, including the fitting of both planar 
surfaces and quadric surfaces to locally approximate the shock wave 
surface. The quadric surface formulation displayed a tendency to 
produce a (local) surface with undulations. The planar surface 
representation did not exhibit this effect, and, for fine mesh spacings, 
produced results essentially the same as the representation that was 
ultimately selected for use in the numerical algorithm. However, 
the accuracy of the planar surface representation suffered at coarse 
mesh spacings. The shock wave surface formulation that was selected 
for use in the algorithm is presented below. 
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The shock wave surface is represented as a family of straight 
lines between two space curves, as illustrated in Figure D.l, The 
space curves represent either the intersection of the shock wave with 
a solution plane (which is a plane of constant x), or the intersec- 
tion of the shock wave with a solid boundary (i.e. , an interplanar 
ring of shock wave solution points). Each space curve is represented 


by the two quadratic 

expressions 



r* (0) = a. 

2 

+ b-9 + c .0 

i i 

(1=1.2) 

(D.2) 

x^e) = d. 

+ e.0 + f.e 2 

(1=1,2) 

(D.3) 


where r. Is the radius of a point on space curve i (i=l,2), x,. is the 
corresponding axial position of a point on space curve i, and 0 is the 
polar angle given by 

0 = tan’^z/y) (0.4) 

where y and z are the coorii nates of a point on the space curve. In equa- 
tions (D.2) and (D.3) t the coefficients a. to f . (i=l,2) are determined 
by fitting these expressions to three known points on each space 
curve as described in Appendix C. When the space curve lies in a 
solution plane, x of course has no 0 dependency. 

Once equations (0.2) and (0.3) are determined for the two space 
curves, the shock wave surface is represented as an infinite family of 
straight lines between the two space curves, where each straight line 
falls in a meridional plane (i.e., a plane of constant 0), Conse- 
quently, for a given value of 0 and x, the shock wave surface is 
represented by the linear interpolation formula 
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FIGURE D. I. SHOCK WAVE SURFACE REPRESENTATION 
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(x - x« ( 9 ) ) (x - X, (0) ) 

r(x - 9) °T » 7 e) ■\( 8 )) r i< e > t - (x 2 Te) i x,(o)) r ;< 9 > < D - 5 > 

In equation (D.5), r(x,0) is the shock wave radius at axial position x 
and polar angle 0, r^(9) and x^ (©) are given by equations (D.2) and 
(0.3), respectively, for one of the space curves, and r 2 (e) and x,,{e) 
are given by equations (D.2) and (D.3), respectively, for the other 
space curve (see Figure D.l). A strong point of this representation 
is that a smooth (local) surface is produced because linear interpola- 
tion is performed for the shock wave radius in a meridional plane, 
while transverse curvature information is introduced through equations 
(D.2) and (D.3). 

I 

4. STREAMLINE- AND BICHARACTERISTIC-SURFACE INTERSECTIONS 

I 

A number of unit processes require determining the intersection 

i 

l 

point of either a streamline with the shock wave, or a bi characteristic 
with either the shock wave or a solid boundary. The technique used 
is the same for all cases and is presented below. 

A streamline or bicharacteristic may be represented by the equa- 
tion 

dx i = T.dt (i=1 ,2,3) (0.6) 

where x,j (1=1 ,2,3) denotes the three cartesian coordinates x, y, and z, 
respectively, and t is a parameter proportional to the length of the 
streamline or bicharacteristic. For a streamline, the parameter r. in 
equation (D.6) is given by 
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0 - 1 . 2 . 3 ) 


(0.7) 


where u^ (1=1 ,2,3) denotes the velocity components u, v, and w, 
respectively. For a bi characteristic, r. is given by 

r i = u. + cre,.cos<{> + c@.jSin4> (i = l,2,3) (0.8) 

where , B. , 4 >, and c are the parameters employed in Butler's 
parameterization of the Mach cone (24), which is discussed in Appendix B. 
Using equation (D.6), the following equation may be written. 

dx/r 1 = dy/r 2 = dz/r 3 (0.9) 

Solving equation (0.9) simultaneously, the linear expressions 


y » [y k - ( r 2 /r,)x k ] + dyr^x 


z ■ t z k ■ ( r 3/ r ])X| ( ] + (Tj/r^x 


( 0 . 10 ) 
( 0 - 11 ) 


may be obtained, where x^, y^, and are the coordinates of a known 
point on the streamline or bicharacteristic, while x, y, and z repre- 


sent the coordinates of the point of intersection of the streamline 
or bicharacteristic with a surface (see Figure D.2). 


An iterative procedure is employed to determine the coordinates 
x, y, and z. First, the values of l\ (1=1 *2,3) are evaluated at the 
known point. Then, a trial value is assumed for the axial coordinate 
x. From equations (D.10) and (D . 1 1 ) , the corresponding coordinates y 
and z may be obtained. Then, the radius r - (y + z ) and the 

_ i 

polar angle 9 = tan (z/y) of the assumed intersection point may be 
computed. From the assumed value for x and the calculated value for e, 
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the body radius r [determined from the tabular wall data or equation 

** 

(D.l)] or the shock wave radius r [given by equation (D.5)] may be 

* *+ 

obtained. The difference between r and r is reduced to within a 
specified tolerance by employing a numerical relaxation technique 
(secant method) which iterates on x. Once convergence has been ob- 
tained, the values of r. at the intersection point are computed using 
the tri variate interpolation method discussed in Appendix G. Appropri- 
ate averages of the values of I\ at the known point and the intersec- 
tion point are then formed, and the entire process is repeated until 
overall convergence is obtained. 

It should be noted that it is possible to use 6, instead of x, 
as the variable upon which the iterative scheme is based. The resulting 
formulation, however, is singular when the streamline or bi character- 
istic lies in a meridional plane. 
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APPENDIX E 


UNIT PROCESSES 


1. INTRODUCTION 

Computation of the flow field requires that a variety of unit 
processes be employed. These subalgorithms may be classified into four 
major types: interior point, solid boundary point, field-shock wave 

point, and body-shock wave point. Computation of the external flow 
field about the forebody portion of the centerbody requires using the 
basic versions of the first three aforementioned algorithms. Computa- 
tion of the internal flow field, with its attendant reflected shock 
wave system, requires using the basic interior point and solid boundary 
point algorithms plus modified versions of these routines, as well as 
the other unit processes. All of the unit processes are presented in 
this appendix. 


2. SUMMARY OF THE CHARACTERISTIC EQUATIONS 

The equations for the characteristic surfaces and the compatibility 
equations valid along these surfaces are developed in Appendix B. A 
summary of the pertinent results is given below. 

For steady three-dimensional supersonic flow, compatibility equa- 
tions may be written which are valid when applied along either stream- 
lines or bi characteristics. A streamline is represented by the equa- 
tion 
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(1-1 ,2,3) 




dx^ = u i dt 


(E.1) 


where (1=1 ,2,3) denotes the three cartesian coordinates x, y, and 
z, respectively, u.. {1 = 1 ,2*3) denotes the corresponding velocity com- 
ponents u, v, and w, respectively, and t is the time of travel of a 

fluid particle along the streamline. The compatibility equations valid 

* 

along a streamline are given by 


dP _ a 2 le. „ 
dt dt 


F 

e 


dP 

it 


+ pu 



i dt 



(E.2) 


(E.3) 


where P denotes the pressure, p is the density, a is the sonic speed, 
(i=l,2,3) denotes the transport forcing terms in the x, y, and z 
component momentum equations, respectively, and F g is the transport 
forcing term in the energy equation. The operator d( )/dt in equa- 
tions (E.2) and (E.3) represents differentiation in the streamline 
direction. The forcing terms F^ and F fi are defined by equations (A. 6) 
and (A. 27), respectively. 

A bi characteristic, which is a ray or generator of the Mach cone, 
is represented by 


dx- = (u i + ca^coso + ce^sinO)dt (i=l,2,3) (E.4) 

where e is a parametric angle denoting a particular element of the Mach 
cone and has the range 0 0 ^ 2n, t is the time of travel of a fluid 


★ 

Repeated Indices imply summation over the range of 1 to 3 unless other- 
wise noted. 
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particle along the streamline that is the axis of the Mach cone, and c 


is defined by 


c 2 . qV/(q 2 - a 2 ) 


(E.5) 


where q is the velocity magnitude. The vectors a. and p. in equation 
(E.4) are parametric unit vectors with a.., 8^, and u^/q (i-1,2,3) 
forming an orthonormal set. The compatibility equation valid along a 
bicharacteristic is given by 


riP *> 

^ + pc(o^cos0 + S^sin©)-^- 3 $ - pe 4 (a.sin© 

9U. 

- 3.cose)(a.sin© - e.cosSk-^- 

« J J oX • 


(E.6) 


In equation (E.6), the operator d( )/dt represents differentiation in 


the bi characteristic direction, and the parameter $ is given by 
» - (c Z /a 2 )(F e - a^F.) 


(E.7) 


where n i is the 1th component of the wave surface unit normal and is 
given by 

o 

n i * (a/c) (cu^/q - axos© - e^sin©} (i=l,2,3) (E.8) 

In addition to the above relations, the following noncharacteristic 
relation is applied along a streamline 


jp p "U* 

SL s 0 . pc (a. a- + 8-6.)v- L 
dt M ' rj p iT3x. 
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where the operator d( }/dt represents differentiation in the streamline 
direction, and the parameter cr is given by 


».attaat' : ■)- . J2Efe££S£ i MS. I 


o = (c 2 /a 2 )F e •■ (c 2 /q 2 )(u j F 1 ) (£.10) 

Equations (E - 1 ) to (E. 10) form the basis of the numerical inte- 
gration method. 

3. GENERAL COMMENTS CONCERNING THE UNIT PROCESSES 

An inverse marching scheme is employed in the numerical algorithm. 
The solution is obtained on space- like planes of constant x, with the 
x-axis being the longitudinal axis of the centerbody and cowl. For the 
internal flow field, the solution is also obtained on the space curves 
which represent the intersection of the internal shock wave with the 
solid boundaries. These space curves are defined by the locus of shock 
wave solution points. 

Except in the vicinity of a shock wave-solid boundary intersec- 
tion, the distance between successive solution planes is determined 
by the application of the Courant-Friedriehs-Lewy (CFL) stability 
criterion, which is presented in Appendix F. The axial step in the 
vicinity of a shock wave-solid boundary intersection is controlled by 
special constraints which are also discussed in Appendix F. 

Each of the unit processes is presented below. In general, a unit 
process is divided into a predictor step and a number of enusing cor- 
rector steps. In most cases, a unit process employs an outer iterative 
loop for determination of the flow properties at the solution point, 
and an inner iterative loop (or loops) for location of bi characteristic- 
initial-value plane intersection points, etc. The terms "inner" and 
"outer" are used in this context in the following discussions. 
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4. INTERIOR POINT UNIT PROCESS 

Figure E.I is a depiction of the computational network used in the 
determination of the solution for a typical interior point. Points (1) 
to (5) are located on the initial-value plane which is a plane of 
constant x on which the solution is known. Points (1) to (4) represent 
the intersection points of four rearward- running bi characteristics with 
the initial-value plane, and point (5) is the intersection point of the 
streamline with this plane. Point (6) is the interior solution point, 
which is located at the intersection of the forward projection of the 
streamline with the solution plane. The axial (x) distance between 
the initial -value plane and the solution plane is determined by either 
the application of the CFL stability criterion, or, in the vicinity 
of a shock wave-solid boundary intersection, by the special constraints 
discussed in Appendix F. 

Interpolated values of the three velocity components u, v, and w, 
the pressure P, and the density p are required at the bi characteristic- 
initial-value plane intersection points, points (1) to (4) in Figure 
E.I. For this purpose, the following bivariate interpolation poly- 
nomial is employed 

f(y,z) » a 1 + a^y + a 3 z + a 4 yz + a 5 y + a & z (E.ll) 

where f(y,z) denotes a general function of the coordinates y and z. 

The coefficients (i°l to 6) in equation (E.ll) are determined by a 
least squares fit of nine data points in the ini tial- value plane 
[point (5) and its eight immediate field point neighbors]. The detailed 
implementation of equation (E.ll) is discussed in Appendix C. 
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FIGURE E.l . INTERIOR POINT COMPUTATIONAL NETWORK 
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In addition to using interpolated values for the flow properties 
at points 0) to {4) in Figure E.l, interpolated values are also employ- 
ed at point (5), the streamline base point, even though this is a field 
solution point. As shown by Ranson, et al. ( 9), this interpolation 
is required to produce a stable numerical scheme. 

The interior point unit process is initiated by locating the 
solution point, point (6). This is accomplished by extending the 
streamline for.ard from point (5) to intersect the solution plane. 

The coordinates of point (6) are obtained using the following finite 
difference form of equation (E.l). 

x f (6) - Xl (5) = ^Cu, (5) + u.(6)][t(6) - t(5)] <i=l,2,3) (E.12) 

In applying equation (E.12) for the predictor {first outer iteration), 
u^{6) is equated to u^(5), whereas, for the corrector {ensuing outer 
iteration), the previously obtained value of u . (6 ) is used. 

Equation (E.12) is first applied for i-1 (i.e., the x-coordinate 
direction). The axial step [x(6) - x{5)] is determined prior to the 
application of the unit process. Hence, the time parameter [t(6) - 
t(5)] may be obtained. Then, equation {E.12) is applied for i=2 and 
i=3 to determine y (6) and z(6). 

At this point, four bi characteristics are extended backward from 
the solution point to intersect the initial-value plane. This is 
accomplished by applying the following finite difference form of 
equation (E.4). 
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x t (6) - x i (k) = ]4[u..(k) + u. (6)] + [c(k) + e(6)][a.Gose(k) 

+ 3 i stne(k)3}J[t(6) - t(k)] (i=l ,2,3) (E.13) 

In equation (E. 1 3) » k denotes the bi characteristic intersection points 
in Figure E.l and has the values 1, 2, 3, and 4 corresponding to the 
e(k) values of 0, tt/2, tt, and 3ir/2, respectively. The bi characteristic 
intersection points are determined in an inner iterative loop. That is, 
for every outer iteration that is performed to determine the flow 
properties at point (6), a number of inner iterations are performed to 
locate points (!) to (4). On the first inner iteration of the predictor 
(the first outer iteration), u- (k) and c(k) are equated to u- (5) and 
c(5), respectively, for each of the four bi characteristics. On ensuing 
inner and outer iterations, the flow properties previously obtained at 
each of the bi characteristic intersection points are used. The flow 
properties at these points are determined by employing the bivariate 
interpolation polynomial given by equation (E.ll). Moreover, as was 
done for equation (E.12), for the predictor (the first outer iteration), 
the flow properties at point (6) in equation (E. 13) are set equal to 
those at point (5), whereas, for the corrector (ensuing outer itera- 
tions), previously computed values of the flow properties are used at 
the solution point. 

Equation (E.13) is first applied for i=l (i.e., the x-coordinate 
direction). The axial step [x(6) - x(k)] is determined prior to the 
application of the unit process. Thus, the time parameter [t(6) - 
t(k)] may be obtained for each of the four bi characteristics. Then, 
equation (E.13) is applied for i=2 and i=3 to determine y(k) and z{k) 
for each bieharacteristic. 
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The parametric unit vectors a. and appearing in equation 
(E.13) are arbitrarily fixed at the solution point, point (6). Butler 
( 24) > in his original work, held and B^ constant along a bichar- 
acteristic but varied 0 In order to insure that the bi characteristic 
remained tangent to the Mach cone. Ransom, et al. (9 ) held 0 constant 
along a bi characteristic but varied and to satisfy this tangency 
condition. As noted by Gline, et al. (25), Butler (41) later realized 
that it is not necessary to satisfy the tangency condition in order 
to achieve second-order accuracy in the resulting overall numerical 
algorithm. As a consequence, in the present analysis, both 0 and the 
unit vectors and $. are held constant along the bicharacteristics. 

For the external flow field integration, ou and B* are selected to 
straddle the projection of the pressure gradient in the initial -value 
plane. For the internal flow field integration, cx^ and B. are chosen 
to straddle the meridional plane. 

Once the positions of and the flow properties at points (1) to (4) 

have been determined for a given outer iteration, the transport forcing 

functions F , F , F , and F are computed at each of these points 
x y 2 G 

and at the streamline base point, point (5), as described in Appendix G. 
Approximations for the transport forcing functions at point (6) are also 
made at this stage as described in Appendix G. The system of non- 
linear compatibility equations is then solved for the flow properties 
at point (6) as outlined below. 

The compatibility equations valid along a streamline are given 
by equations ( E . 2) and (E.3). Writing those relations in finite 
difference form yields 
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(E.14) 


CP(6) - P(5)]/[t(6) - t(5) ] - |[a 2 (5) + a 2 (6)][p(6) 
- p(5)]/[t(6) - t(5)] = i{F e (5) * F e (6)] 


[P<6) - P(5)]/[t(6) - t(5)] + ^Cp(5)u i (5) + p(6) Ui (6)][u.(6) 

- u,(5)]/[t(6) - t(S)] * |Cu 1 {E)F 1 <5) + u i (6)F 1 (6)] (E.15) 

The noncharacteristic equation, given by equation {E.9)» is also 

applied along a streamline. Writing that equation in finite differ- 
ence form gives 

[P(6) - P(5}]/[t{6) - t(5)] = |Ca(5) + a( 6)] 

- ■|p(5)c 2 (5)(a i a ; . + e^^)3u^/3Xj (5) 

- •^(6)c 2 {6)(a i a J . + (E. 16) 

In equation (E.16), o is given by equation (E.10), and 9u./9x.(k) de- 

* J 

notes the appropriate partial derivative evaluated at point ( k) in 
Figure E.l. Partial derivatives taken with respect to y and z are 
found by analytically differentiating equation (E . 1 1 ) . Partial deriva- 
tives taken with respect to x are then found by using the governing 
partial differential equations. 

The compatibility equation valid along a bi characteristic is given 
by equation (E.6). For 0 values of 0, -rr/2, tt, and 3tr/2, equation (E.6) 
becomes 


dP . „ 
(jtj' + p co < 


du. 2 9u. 

i dt^ = $ 1 * pe 6 i^j 9x7 


(£. 17 ) 
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dP 


du. 


9u. 


ir- + peg. tt- = - pc a. a , 

dto i at, 2 K i j 9x- 


(E, 18) 


dP du 1 2 3u i 

dt 3 ' pca 1 dt 3 ' *3 ‘ pC e i B j 3Xj 


(E.19) 


dP du i 2 3u i 

dt^ • pc6 i dt^“ *4 ' pc Vj 5xJ 


(E.20) 


In equations (E.17) to (E.20), the operator d( }/dt k denotes differenti- 
ation along the bi characteristic corresponding to © (k ) , and is de- 
termined from equation (E.7). Writing equations (E.17) to (E.20) in 


finite difference form yields 


[P(6) - P(l)]/[t(6) - t(l )] + ±te(l)c(l) 

+ p(6)c(6 )]ot i [u 1 (6) - u, (l)]/[t(6) - t(l)] 

- |C^(1) + (6)1 - |p(l)e 2 (l)e.3 j 3u i /9x j (l) 

- ^p(6)c 2 (6)@ i 6- i 3u i /3x j .(6) (E.21) 

[P(6) - P(2)]/[t{6) - t(2)] + “Cp(2)c(2) 

+ p(6)e(6)36 i llu i (6> - u.(2)]/[t(6) - t(2)] 

= ^C* 2 (2) + « 2 (6)] - |p(2)c 2 (2)a i a j 9u./3x j (2) 

- ^p(6)c 2 (6)a^cXj3u^/3Xj(6) (E.22) 
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[P(6) - P(3)]/[t(6) - t(3)] -ltp(3)c(3) 

+ p( 6 >c (6)301^ [u^ ( 6 ) - u i (3)l/[t(6) - t(3)] 

■ |t* 3 (3) + * 3 ( 6 )] - ip(3)c Z (3)S i 6 j »u i /ax J <3) 

- ip(6>e 2 (6)e i 6 j ati i /3x j (6) (E.23) 

[P(6) - P(4)]/[t(6) - t(4)] - ±£p(4)c(4) 

+ i»(6)c(6>]B 1 Cu 1 (6) - u.(4)]/[t<6) - t(4)] 

• g[* 4 (4) + 4 4 (6)] - 5P(4)c 2 (4)o.oi ;i 3u./3x j (4) 

- |p(6)c Z (6)a 1 a j 3u./3x j (6) (E.24) 

It was noted in Appendix B that only three wave surface compatibil 
ity relations are independent. To obtain three independent relations, 
linear combinations of equations (E.21) to (E.24) and the nonchar- 
acteristic relation, equation (E.16), are formed in such a manner as 
to algebraically eliminate the cross derivative terms at the solution 
point [i.e., terms containing 3u./3x.(6)]. Subtracting equation (E.23) 

* iJ 

from equation (E.21) yields 
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t+ywr . -tcarrart ■sxrtt&fsn * tosurW 'JftJlKrr-t ’PWskjti- 


[P(6) - P(1)3/[t(6) - t(l )2 - [P(6) - P(3)]/[t(6) - t(3)] 

+ j{p(1)c( 1) + p(6)c(6)]o.£u j (6) - u.(l)]/[t(6) - t(l)3 

+ |tp(3)c(3) + p(6)e(6)]a^[u^(6) - u f (3)Mt(6) - t(3)] 

» |£*,(i) + t,(6)3 - |f> 3 (3) + 4. 3 (6)] 

- ^Ote^n^BjSu^XjO) + ip(3)c 2 (3)B 1 B j 9u 1 /8x j .(3) (E.2S) 

Subtracting equation (E.24) from equation (E.22) yields 

[P(6> - P(2))/[t(6) - t(2)] - [P(6) - P(4)]/[t(6) - t(4)] 

+ |{p(2)c(2) + 0 ( 6 ) 0 ( 6 ) 36 ^. (6) - u.(2)3/[t(6) - t(2)3 

+ |(o(4)c(4) + p(6)c(6)36 1 -Cu i (6) - u. (4)3/[t(6) - t(4)3 

- ^£* 2 ( 2 ) + + *q(6)3 

- ^o(2)o 2 (2)a ( oij3u./3Xj(2) + |p(4)c 2 (4)a 1 -aj9u j /3x ; .(4) (E.26) 

Adding equations (E. 21 ) and (E.22) and subtracting equation (E.16) from 
the sum yields 

[P(6) - P{l)]/[t(6) - t(l)] + [P(6) - P(2)]/[t(6) - t(2)] 

- [P{6) - P(5)]/[t(6) - t(S)] 

+ |(o(l)c(l) + p(6)c(6)]cx i Cu 1 (6> - U.(l)3/[t(6) - t(l)] 

+ i(p(2)c(2) * pfejcfenSjU^) - Ui (2)3/[t(6) - t(2)3 
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- ^c* 1 (l) + ^(6)] + |C$ 2 (2) + $ 2 (6)] - lfcr(5) + o<6)] 

- |p(l)c 2 (l)e i e j au i /9x j (l) - |p(2)c 2 (2)' J i i a j 9u 1 /3x j (2) 

+ |p(5)e 2 (5)(Gt i a^ + g^g.)3u^/9Xj(5) ( E . 27 ) 

Equations (E.14), (E.1S), (E.25), (E.26), and (E.27) are the five 
finite difference equations which are used to solve for the flow 
properties u(6), v (6 } , w(6), P(6), and p(6). Since these equations are 
nonlinear, an iterative scheme is required to obtain the solution. On 
the first outer iteration (the predictor), all of the flow properties 
at point (6) appearing in the coefficients of the derivatives in the 
above set of equations are set equal to the respective properties at 
point (5). This produces a system of simultaneous linear equations 
which is solved using a Gaussian elimination method with complete 
pivoting (40). On ensuing corrector applications (outer iterations), 
previously computed values for the flow properties at point (6) are 
employed in the scheme. This method is similar to the Euler predictor- 
corrector algorithm used to obtain the solution for initial -value 
problems for ordinary differential equations, and can be shown to have 
second-order accuracy either by direct numerical calculation ( 9 ) or 
by substituting an exact solution into the difference equations and 
expanding the resulting terms in a Taylor series and thereby determining 
the truncation error. The iterative scheme is terminated when all five 
flow properties at point (6) have converged to within specified toler- 
ances. 


as* 


161 


5. SOLID BOUNDARY POINT UNIT PROCESS 

Figure E.2 is a depiction of the computational network used in 
determining the solution for a typical point on a solid boundary. The 
point notation used in Figure E.2 is the same as that used in Figure 
E.l (interior point scheme). In this unit process, however, point (4), 
corresponding to the bi characteristic with 0 - 3 tt/ 2, falls outside of 
the flow field and cannot be employed. Furthermore, the streamline 
points (5) and (6) lie on the stream surface formed by the solid boun- 
dary. The formulations used for representing the solid boundaries 
are presented in Appendix D, 

The boundary condition used in this unit process is simply that 
the flow be tangent to the surface of the boundary at the solution 
point, point (6) in Figure E.2. Let n b * (1*1 ,2,3) denote the x, y, 
and z components, respectively, of the outward unit normal to the solid 
boundary surface. Then, the flow tangency boundary condition may be 
written as 

u i (6) n bi (6) = 0 (E.28) 

The solid boundary point unit process is virtually identical to 
the interior point unit process, except that the wave surface compati- 
bility equation valid along the bi characteristic corresponding to 0 = 
3 tt/ 2 is not employed. That equation is replaced by equation (E.28). 
Thus, the system of compatibility equations used for determining the 
solution at a solid boundary point consists of equations {E. 14) , (E.15), 
(E.25), (E.27), and (E.28). This system of equations is solved using 
the same iterative scheme that was employed in the interior point 
solution. 
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The location of the solution point, point (6) in Figure E.2, ob- 
tained by applying the finite difference form of the streamline equa- 
tion, equation ( E. 1 2) , is adjusted along the projection of the body 
normal in the solution plane so that the solution point lies on the 
solid boundary. The orientation of the parametric unit vectors and 
is selected such that g. * - n fai (i=l,2,3), and ou (i=l,2,3) is found 
by employing the orthonormal relations between oc., 3^ , and u./q. This 
selection for the reference vector set produces a computational network 
in which the bi characteristics corresponding to © = 0, tt/2, and tt 
intersect the initial -value plane for convex boundaries. For concave 
boundaries, those bi characteristics intersect an extrapolation of the 
initial -value plane (the required extrapolation is assumed to have an 
error third-order in step size). The bi characteristics corresponding 
to 6 - 0 and tt lie in the elemental plane which is tangent to the solid 
boundary at point (6). 

6. BOW SHOCK WAVE POINT UNIT PROCESS 

A depiction of the computational network used in determining the 
solution for a typical bow shock wave point is given in Figure E.3. A 
segment of the shock wave surface extending from the initial-value plane 
to the solution plane is shown in this figure. The space curve (A) is 
defined by the intersection of the shock wave with the ini tial- value 
plane, whereas, space curve (B) is defined by the intersection of the 
shock wave with the solution plane. The axial distance between the 
ini tial -value plane and the solution plane is determined by the appli- 
cation of the CFL stability criterion. 

The bow shock wave solution point is denoted by point (2) in 
Figure E.3. The flow properties upstream of the shock wave are known 

atMtumamasniM ' 







a priori . Hence, In the following discussion, the flow properties u{2), 
v(2), w(2), P{2), and p (2) refer to the properties at point (2) down- 
stream of the shock wave. Point (1) is the intersection point of a 
rearward- running bicharacteristic with the initial -value plane. This 
bi characteristic is extended backward from the solution point. Point 

(3) is an interior point in the solution plane which is used to define 
the meridional plane in which the shock wave solution point lies. Point 

(4) is the intersection point of space curve (A) with the meridional 
plane which passes through points (2) and (3). 

In this unit process, a local cartesian coordinate system is 
employed for the description of the orientation of the local shock wave 
surface. This local coordinate system has coordinates x', y'» and z‘, 
where x' is coincident with the x-axis, y* is in the radial direction 
corresponding to the meridional plane which subtends an angle 6 with the 
fx,y)-plane» and z* is normal to the (x* ,y‘ )-plane (see Figure E.3). The 

A A A 

unit vectors in the x, y, and z directions are denoted by i, j, and k, 
respectively, whereas, the unit vectors in the x', y'» and z' directions 

A A A 

are denoted by i', j', and k', respectively. A vector quantity A may be 
represented in these coordinate systems by 


A = A i + A j + A k (E . 29) 

A j 4 

J = A x ,i‘ + Ay, j * + A z ,k' (E . 30} 

The relationships between the respective components in equations (E.29) 
and (E.30) are given by 



{E.31J 
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• .V.v*r«-. ,- '.!>■-•• <: • .: .-,*... :- 




A , = A cose + A sine 
y j 2 


A z , = A z cose - Arsine 


\ ■ V 


A = A .cose - A , sine 

J J “ 


A z * A z ,cose + Arsine 


{ E - 32) 
(E. 33) 
(E.34) 
(E.35) 
(E.36) 


The orientation of the local shock wave surface is specified by 
a set of unit vectors referenced to the (x‘ ,y 1 »z' )-system. This set of 

A 

unit vectors, illustrated in Figure E.4, consists of a unit vector n $ 

A A 

which is normal to the shock wave surface and two unit vectors JL and t 

A 

which are tangent to this surface. The tangential unit vector t lies 
in the meridional plane [(x 1 ,y* )-plane], subtends an angle (j> with the 
x'-axis, and is defined by the intersection of the shock wave with the 

A 

meridional plane at point (P). The tangential unit vector j, lies in 
the transverse plane [(y ' ,z' )-plane], subtends an angle a with the 
z'-axis, and is defined by the intersection of the shock wave with the 

A A 

transverse plane at point (P). The tangential vectors t and i are 
therefore given by 


t = cosep i ' + sin<f> j’ 


(E- 37) 


JL = sina j' + cosa k' 


(E.38) 


The shock wave normal unit vector, denoted by n g , is given by 

AAA A 

n $ = JL x t/|JL x t| 


(E.39) 
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FIGURE E.4. UNIT VECTORS FOR SPECIFICATION OF SHOCK WAVE 

SURFACE ORIENTATION 



v :-:| MII'fTl ., - .... , %v ^ - if --r{ ...r I. i.l . li ■ I . 

I The interior point and solid boundary point unit processes achieve 

second-order accuracy by using local iteration. In local iteration, 
a corrector application employs previously determined flow property 
values at the solution point, but does not require using flow property 
values at other points in the solution plane. The shock wave point 
unit process, however, requires that global iteration be performed in 
order to achieve second-order accuracy. In global iteration, a cor- 
rector application employs previously determined flow property values 
not only at the solution point, but also at neighboring points in the 
solution plane. As a consequence, before a corrector application in 
global iteration can be performed, the entire solution plane {or at 
least an appropriate section of it) must be determined by a prior 
calculation. In practice, since the interior point and solid boundary 
point schemes require local iteration only, the interior point and 
solid boundary points are computed first. Then, a prediction for each 

shock wave solution point is made, thereby giving a tentative solution 

for all of the shock wave points. Then, a global iteration is con- 
ducted for the shock wave solution points using the previously de- 

termined field points in the solution plane. In the following discus- 
sion, the term "predictor” will refer to the first application of the 
shock wave point unit process used to obtain an initial estimate of the 
solution without using field point data in the solution plane. The 
term "global corrector" will refer to the application of the shock wave 
point unit process which uses field point data in the solution plane. 
The shock wave point unit process is now outlined. 


The shock wave point unit process is initiated by locating the 
solution point, point (2) in Figure E.3. The meridional plane in which 
the solution point lies is arbitrarily selected to contain point (3). 
Point (3) is the interior solution point adjacent to the shock wave sur- 
face whose location is determined prior to the application of the shock 
wave point unit process. The angle subtended by a meridional plane 
and the (x,y)-plane is denoted by 0. Then 

6(2) > 6(3) * tarf 1 [z(3)/y(3)] (E.40) 

Denote the radial position of a point by r. Then the radial position 
of point (2) Is obtained from 

r(2) = r(4) + [x(2) - x(4)] tan |{f>(2) + 4(4)]} (E.41) 

where [x(2) - x( 4)3 is the axial distance between the initial-value 
plane and the solution plane and is determined by the GFL stability 
criterion. On the first application of equation (E.41), the shock wave 
angle $(2) is equated to 4>(4), whereas, on ensuing applications, the 

previously determined value of 4>(2) is used. At point (4), the radial 

position r{4) and shock wave angle <f>(4) are determined by interpolation 
using the quadratic univariate formulae 

r(e) = a 1 + a 2 e + a 3 e 2 (E. 42) 

<J>(0) = b ] + b 2 © + b 3 0 2 (E.43) 

In equations (E.42) and (E.43), the coefficients a. (i=l,2»3) and b^ 
(i=l,2,3) are determined by fitting these expressions to three local 
shock wave solution points on space curve (A) as described in Appendix C. 


For the ease of axi symmetric flow, or on a plane of flow symmetry in 
three-dimensional flow, point (4) coincides with a previously determined 
shock wave solution point so the interpolation would not be required. 

In general, however, point (4) does not coincide with a known point 
so the interpolation is necessary. 

After the solution point has been located, the shock wave normal 
unit vector n s at the solution point is found by forming the normalized 

A A 

cross product of the tangential unit vectors l and t [see equation 

Aw 

( E . 39 ) 3 . The tangential vector t is obtained by usiij the current 

A 

value of (2 ) in equation (E.37). The tangential vector a is obtained 
by using the current value of a(2) in equation (E.38). For either 
space curve (A) or space curve (B), the value of a(2) may be obtained 
from 


<*(2) = tan' 


1 fl_ dr 
r d0 


0 ( 2 ) 


(E.44) 


For a predictor application, the analytical form of r(e) used in equa- 
tion (£.44) is given by equation (E.42) applied along space curve (A), 
whereas, for a global corrector application, r(©) is obtained from 
equation (E.42) applied along space curve (B). 

After the shock wave normal unit vector has been determined, the 
local Hugoniot equations may be applied across the shock wave, thereby 
yielding a solution for the flow properties u(2), v(2), w(2), P (2 ) , and 
p{2). In general, the local Hugoniot equations take the form (5) 


p u V nu = °d V nd 

P + p V 2 ., • P. + p.V 2 , 
u u nu d H d nd 


(E.45) 
( E . 46 ) 
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(E.47) 


*»■ * v id 

(E.48) 

h u + - h d + «d /2 

(E.49) 

h - h(P,p) 

(E.50) 


In equations (E.45) to (E.50), h is the enthalpy per unit mass, q is the 
velocity magnitude (q = u + v + w },V n is the velocity component in 

A A *“ 

the ~n s direction, is the velocity component in the t direction, 

A 

is the velocity component in the i direction, and the subscripts u and 
d denote the properties on the upstream and downstream sides of the 
shock wave, repsectively. Equations (E.45) to (E.50) are solved simul- 
taneously for the downstream flow properties. To obtain the velocity 
components V nu » , and V £u , the upstream velocity vector is first 
transformed from the (x>y,z)-system to the (x * »y‘ ,z' )-system using 
equations (E.31) to (E . 33) , after which the appropriate dot products 

A A A 

are formed with -n s » t, and z. Similarly, the downstream velocity 
components V^, V t(j , and V ^ are transformed back to the (x,y,z)-system 
after the local Hugoniot equations have been applied. 

In the computer program, the local Hugoniot equations are contained 
in a separate subroutine. The assumed thermodynamic model is that of 
a thermally and calorically perfect gas. Other thermodynamic models 
may be used by suitably modifying the existing subroutine or replacing 
it. For the assumed model of a thermally and calorically perfect gas, 
the pressure ratio across the shock wave is given by 


172 




, 2y m 2 . Y - 1 
P u y + 1 nu y + 1 

where M nu is the incident normal Mach number given by 




/a 


u 


(£. 51 ) 


(E.52) 


and y is the specific heat ratio. Using the result of equation (E.51), 
the density ratio across the shock wave is given by 


P d (y + 1)/(Y - 1) + (P u /P d ) 

P u = i '+ L(y + D/(y - i)](P u /P d ) 


With the downstream pressure and density determined, the downstream 
normal velocity component V ncl may be obtained from equation (E.46), and 
the tangential downstream velocity components and may be com- 
puted from equations (E.47) and (E.48). Transformation of the down- 
stream velocity components back into the {x,y,z)-system yields the 
required flow properties at the solution point. 

At this stage, a rearward-running bi characteristic is extended from 
the solution point, point (2), back to the initial-value plane, inter- 
secting this plane at point (1), as illustrated in Figure E.3. This 
is accomplished by employing the following finite difference form of 
equation (E.4) evaluated for the parametric angle 9 » tt/2. 


(,(2) - x.(l) - | {[11,(1) + u, (2)] 


+ [c(l ) + c(Z)Js. 


kt(2) - ton (i-i. z, 3) 


(E.54) 


As in the interior point and solid boundary point schemes, an inner 
iteration is performed to locate point (1). On the first application 
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of equation (E.54), the flow properties at point (1) are equated to 
those at point (2), whereas, on ensuing applications, previously ob- 
tained values of the flow properties at point (1) are used. The flow 
property values at point (1) are found by employing the bivariate inter- 
polation polynomial given by equation (E.ll). The coefficients in 
equation (E.ll) are obtained by a least squares fit of nine data 
points in the initial-value plane using a boundary- type stencil as 
described in Append i C. 

Equation (E.54) is first applied for i=l (i.e., the x-coordinate 
direction). Since the axial step [x(2) - x(l)] is known from the 
application of the CFL stability criterion, the time parameter 
[t(2) - t(l)] may be determined. Then, equation (E.54) is applied for 
i-2 and i-3 to determine y(l) and z(l). For axi symmeti rc flow, or for 
a plane of flow symmetry in three-dimensional flow, point (f) lies in 
the meridional plane which contains points (2) and (3). In general, 
however, for other flow situations, point (1) lies outside of this 
plane. 

The orientation of the parametric unit vector 8^ in equation (E.54) 
is arbitrarily selected such that 

& 3 /e 2 - tan[©(2)] (E.55) 

This relation, in conjunction with the orthonormality conditions 

8^.(2) = 0 (E.56) 

= 1 ( E . 57 ) 

allows the values of 8^ (1=1 ,2,3) to be determined. Since equation 
174 








(E.57) is a quadratic equation, a multiplicity of roots exist for the 
(i=l,2,3). The roots are chosen such that point (1) lies under- 
neath the shock wave in the initial-Vilue plane. Once the values of 
(i-1,2, 3) are determined, the values of (i=l,2,3) are found 
through use of the orthogonality relation between a.. , 3- , and u^/q 

A A 

(i .e. , a = 0 x V/q) . 

After the position of and the flow properties at point (1) have 

been determined, the transport forcing functions F , F , F , and F are 

x y z e 

computed at point (1) as described in Appendix G. Approximations for 
the transport forcing functions are also made at point (2) at this 
time as described in Appendix G. 

At this stage, the wave surface compatibility equation correspond- 
ing to the parametric angle s3 = v/2 is applied between points (1) and 
(2). From equation (£.6), the appropriate equation is 


If + ocs < 


du 


3u, 


i dt ^n/2 pC *i l j ax . 


(6.58) 


where $ r j2 stained from equation (E.7) for the parametric angle 6 = 
tt/2 . Writing equation (E.58) in finite difference form, solving for 
the pressure at point (2), and denoting this pressure by P*(2), the 
following equation is obtained. 


P*(2) * P(l) + yO^O) + $ TT/2 (2)]:[t{2) - t(l)] 

- |[p(l)c 2 (l)a i a j 3u./3x j (l) 

+ p{2)c 2 (2)a i a j 3u i /3x j (2)][t(2) - t(l)] 

- {(p(l)c(l) + p(2)c(2)]0.[u i {2) - u.(l)] (E.59) 
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Note that the cross-derivative terms [3u./9x.(k)] in equation 

* J 

(E.59) appear at both point (1) in the initial -value plane and at 
point (2) in the solution plane. In general, these terms can be 
evaluated by employing equation (E.ll) fit to nine data points in the 
appropriate plane, differentiating this expression analytically to ob- 
tain partial derivatives with respect to y and z, and then using the 
governing partial differential equations to obtain the required partial 
derivatives with respect to x. On the predictor application of the 
shock wave point unit process, the flow property field in the solution 
plane is not known, so the cross-derivatives at point (2) are set equal 
to those at point {!). On a global corrector application of the 
shock wave point unit process, the cross deriatives at point (2) are 
evaluated in the manner just described. 

The pressure P(2) is calculated from the local Hugoniot equations. 
The pressure P*(2) is calculated from equation (E.59). The difference 
between P(2) and P*(2) is driven to within a specified tolerance of 
zero by employing a one-dimensional secant iteration scheme which 
iterates on the shock wave angle d> (2 ) . Two initial estimates of 
are required to initiate the subiteration. 

The shock wave point unit process is first applied as a predictor 
for each shock wave solution point. In this application, the value of 
a used in equation (E.38) is obtained by curve fitting points along 
space curve (A), and the cross-derivative terms at the shock wave solu- 
tion point are equated to those terms at the bicharacteristic base 
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point in the initial-value plane, point (1). After a tentative solution 
is obtained for all of the shock wave points, a number of global cor- 
rector applications are performed. Here, the value of a used in equa- 
tion (E.38) is based on data along space curve (B), and the cross- 
derivative terms at the shock wave solution point are evaluated at that 
point. The resulting overall scheme has second-order accuracy when 
the global correction is performed. The global iteration is terminated 
when successive values of » have converged at each of the shock wave 
solution points . 

In the course of the program development, an alternative algorithm 
to the one just presented was devised in an attempt to compute the bow 
shock wave solution points. In this alternative scheme, a multiplicity 
of bicharacteristics were used, and, like the interior point or solid 
boundary point unit processes, linear combinations of the wave surface 
compatibility equations were formed as to algebraically eliminate the 
cross -deri vati ve terms at the solution point. A two-dimensional 
Newton- Raphson method was devised for determining the angles $ and 
explicitly, and second-order accuracy was achieved without resorting to 
global correction. This scheme was successful in computing axisymmetric 
flows, but an apparent instability arose when attempting to compute 
three-dimensional flow fields. 

7. SOLID BODY-SHOCK WAVE POINT UNIT PROCESS 

The solid body-shock wave point unit process is used’ to determine 
the flow properties downstream of the shock wave at a point where the 
shock wave intersects a solid boundary. This unit process is used to 
determine the solution for the points on the cowl on the downstream side 


177 


of the cowl lip shock wave, and for the points on the centerbody or cowl 
on the downstream side of an internal reflected shock wave. The method 
of computation is essentially the same for either application and is 
discussed below. The solution points on the downstream side of the 
incident shock wave at an internal shock wave reflection are computed 
using the field-shock wave point unit process which is presented later. 

A depiction of the computational network used in the solid body- 
shock wave point unit process is presented in Figure E.5. A typical 
solid body- shock wave solution point is denoted by point (P) in this 
figure. At point (P ) , the outward unit normal vector to the solid 
boundary is denoted by n^. The locus of solid body-shock wave solution 
points represents the intersection of the shock wave with the solid 
boundary, and defines space curve (A) in Figure E.5. The intersection 
of the shock wave with the meridional plane passing through point (P) 
is denoted by space curve (B). The tangential unit vectors to space 
curves (A) and (B) at point (P) are denoted by l and t, respectively. 

The unit normal vector to the shock wave at point (P) is denoted by 

V 

As was done for the bow shock wave point unit process, the init 

vectors t, and n are referenced to a local cartesian coordinate 

s 

system {x',y',;’'), where again x' is coincident with the x-axis, y' is 
in the radial direction along the meridian which subtends the angle »* 
with the {x,y)-plane, and z‘ is normal to the (x 1 ,y ' )-plane. The rela- 
tions between the components of a vector in the (x,y ,z }- system and in 
the (x 1 ,y' ,z’ )-systen are given by equations (E.31) to (E.36). As in 
the bow shock wave point unit process, the tangential unit vector t lies 
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MERIDIONAL PLANE 



FIGURE E.5. SOLID BODY-SHOCK WAVE POINT COMPUTATIONAL 

NETWORK 



I 


in the meridional plane C(x’ ,y* )-plane] and subtends the angle 41 with 
the x'-axis. Hence, 


t = costf) i ' + sirxj) j ' 


(E.60) 


Unlike the bow shock wave point unit process, however, the tangential 

A 

unit vector % does not, in general, lie in the transverse plane 
C (y ' tZ* )- plane], but rather it may have a nonzero x' -component. This 
tangential vector along space curve (A) may be represented by 


* ds ds J ds 


(E.61) 


where ds is the differential arc length given by 


(ds) 2 = (dx') 2 + (dy * ) 2 + (dz’) 2 


(E. 62) 


The derivatives in equation (E.61) are obtained by analytically 
differentiating the expressions 


x 1 (© ) = a.j + a^e + a 3 £ 


(E - 63) 


y' (0) = b 1 + b 2 e + b 3 e' 


(E.64) 


z * (© ) = c^ + c^e + c-f 


(E.65) 


In equations (E.63) to (E.65), the coefficients a., b.., and c. 

(1=1 ,2,3) ere obtained by fitting the respective expressions to three 
points on space curve (A) as described in Appendix G. For the cowl lip 
shock wave poi,,cs, space curve (A) is defined by the cowl lip itself 
since the shock wave is assumed to be attached to the cowl lip. In 
this case, the x'-component in equation (E.61) is identically zero. 
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and, as a consequence, £ lies in the transverse plane. Furthermore, if 
the cowl is axl symmetric, the y* -component is also identically zero. 
Alternatively, for computing the downstream properties at a reflected 
internal shock wave, space curve (A) is defined by the intersection 
of the incident shock wave with the solid boundary. Except for an axi- 
symmetric flow field, or for a point on a plane of flow symnetry in 
three-dimensional flow, the x'-component in equation { E . 61 ) is nonzero. 
With the tangential unit vectors determined, the shock wave normal unit 

a 

vector n is obtained from equation (E.39). 

The solid body-shock wave point unit process is initiated by 

A 

determining the body normal unit vector n b and the tangential unit 

A 

vector z at point (P), expressing both of these vectors in the 
(x' ,y ' ,z' )-system. Then, an initial estimate is made for the value of 
$ in equation (E.60), and, by use of equation (E.39), the shock wave 
normal unit vector is obtained. In exactly the same manner as was done 
in the bow shock wave point unit process, the downstream flow properties 
at point (P) are computed by use of equa ons (E.45) to (E.53). At 
this stage, the velocity normal to the body V nb at point (P) is computed 
from the equation 


nb = u dV + V dV * Vbz 


where u^, v^, and w^ are the downstream velocity components at point 
(P), and n bx ,, , , and n b2 , are the components of the body normal unit 
vector, both vectors being expressed in terms of the (x',y',z') 


coordinates. The body normal velocity V nb is reduced to within a 
specified tolerance of zero by varying the angle 4> using a ene- 
dimensional secant iteration procedure. Two initial estimates of 4> are 


required for starting the iterative procedure. Once convergence has 
been obtained, the downstream velocity components are transformed back 
into the (x,y, z) -coordinates using equations (E.34) to (E.36). 

In the course of the program development, an alternative algorithm 
to the one just presented was devised to compute the solid body-shock 
wave points. That algorithm determined the shock normal vector (and 
thereby the downstream properties) by employing the shock wave rela- 
tions which link the flow turning angle and the shock wave angle, both 
these angles being measured from the approach streamline direction in 
a plane defined by the approach velocity vector and the shock wave 
normal vector. Since the shock wave normal vector is required to de- 
fine this plane, an iterative procedure for determining that vector is 
required in this method. This method was tested and produced results 
identical to the method described earlier. However, due to the greater 
complexity of the alternate method, it was not selected for use in the 
final algorithm. 

8. SHOCK-MODIFIED INTERIOR POINT UNIT PROCESSES 

In some situations during the computation of the internal flow, 
the interior point unit process must be applied in a modified form. 

One such application is illustrated in Figure E.6. In this situation, 
the Mach cone, with apex at the solution point, intersects not only the 
initial-value plane but also a solid boundary and an internal shock 
wave. The point notation used in Figure E.6 is the same as that used 
in the computational network of the basic interior point scheme, which 
is illustrated in Figure E.l. The solution point, denoted by point (6) 
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FIGURE E.6. SHOCK- MODIFIED INTERIOR POINT COMPUTATIONAL 
NETWORK (STREAMLINE BASE POINT ON 
INITIAL -VALUE PLANE) 




in Figure E. 6, lies on the current solution plane. Point (5) represents 
the streamline base point on the initial -value plane. As in the basic 
interior point unit process, points (1) to (4) represent the bichar- 
acteristic base points. Point (1), in this case, lies on the surface 
of the internal shock wave, and point (3) lies on the solid boundary. 
Points (2) and (4) lie on the initial-value plane. 

The axial distance between the initial -value plane and the solu- 
tion plane is determined by either the CFL stability criterion or by 
the special constraints which apply when an internal shock wave inter- 
sects a solid boundary. Those procedures are discussed in Appendix F. 

In either case, the axial step is determined prior to the application 
of the unit processes. 

In the overall algorithm for the computation of the internal flow, 
the order of integration is selected so that the shock wave solution 
points and the body solution points are determined before any attempt 
is made to obtain the solution at any of the interior field points 
which lie in the flow field sector that is downstream of the shock wave. 
As a consequence, the flow property fields on the downstream side of 
the shock wave and on the stream surface formed by the solid boundary 
are determined before the solution at an interior point, such as 
point (6) in Figure E.6, is attempted. 

The procedure used to obtain the solution at point (6) in Figure 
E.6 is almost identical to the basic interior point unit process, which 
is presented in Section 4 of this appendix. The major difference be- 
tween the two algorithms is that, in the present case, the bi character- 
istic intersection points on the shock wave [point (1)] and on the 
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solid boundary [point (3)] must be determined in addition to those bi- 
characteristic intersection points [points (2) and (4)] on the initial- 
value plane. Along with the location of these points, flow property 
values and first partial derivatives of the flow properties at these 
points must also be obtained. 

As in the basic interior point unit process, flow property values 
at points (2) and (4) on the initial-value plane are obtained using 
the bivariate interpolation polynomial given by equation (E.ll). The 
coefficients in this equation are determined by a least squares fit of 
nine data points in the initial-value plane as discussed in Appendix C. 
Flow property values at point (1) on the shock wave surface or at 
point (3) on the solid boundary surface are obtained using the tri- 
variate interpolation polynomial 

2 2 

f(x,y,z) = a 1 + a^y + a 3 z + a 4 yz + a 5 y + a g z 

+ a^xy + agXz (E.67) 

The coefficients (1=1 to 8) in equation (E.67) are determined by a 
least squares fit of fourteen data points on either the downstream side 
of the shock wave for interpolation on that surface, or on the solid 
boundary for interpolation on that surface. The detailed implementa- 
tion of equation (E.67) is presented in Appendix C. 

An outline of the unit process used to determine the solution at 
point (6) in Figure E.6 is now presented. The computation is initiated 
by determining the location of the solution point, point (6), using 
equation (E.12) in a manner identical to the procedure employed In the 
basic interior point unit process. After the position of the solution 
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point has been obtained for a given outer iteration, the four bi- 
characteri sties, corresponding to the values of the parametric angle 
© = 0, ti/2, tt, and 3 tt/ 2 in equation (E.13), are extended rearward from 
the solution point to the initial -value plane. From the bicharacteris- 
tic-initial-value plane intersection point coordinates, denoted by 
y*(k) and z*(k) (k=l to 4), the radius r*(k) « [y*(k) 2 + z*(k) 2 ]^ 2 and 
the polar angle 0*(k) = tan 1 [z*(k)/y*(k)] of each intersection point 
are computed. The radius r*(k) is then compared to the shock wave 
radius r & and the body radius r b in the meridional plane defined by the 
polar angle 0*{k). The shock wave radius is determined from the uni- 
variate interpolation polynomial 

r s (e ) = a-j + a^Q + a^e 2 (E.68) 

where the coefficients a. (i=l,2,3) are determined by fitting this ex- 
pression to three shock wave solution points in the initial -value plane 
as described in Appendix C. The solid body radius r b is obtained by 
employing the formulations presented in Appendix D. For the orientation 
shown in Figure E.6, if r & < r*(k) < r^, the bicharacten:, tic inter- 
sects the initial-value plane and the analysis proceeds as in the basic 
interior point unit process. If r*(k) < r , the bicharacteristic 
intersects the internal shock wave. In this ease, the bicharacteristic 
base point location on the surface of the shock wave is found by 
employing the bicharacteristic-surface intersection scheme presented in 
Appendix D. For a shock wave intersection, that scheme requires that 
equation (E.68) also be fitted to three shock wave solution points in 
the current solution plane. If r*(k) > r^, the bicharacten s tic 
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intersects the solid boundary. The bi characteristic base point loca- 
tion on the solid boundary is also obtained by using the iterative 
scheme presented in Appendix 0 As in the basic interior point unit 
process, an inner iteration is performed for locating points (1) to (4). 
Interpolated values of the flow properties at the respective points 
are obtained by using either equation (£.11) or equation (E. 67) * 
whichever is applicable. 

After the bi characteristic base points, points (1) to (4), have 
been located, the first partial derivatives of the flow properties 
with respect to y and z at these points are obtained by analytically 
differentiating the appropriate interpolation polynomial. In a like 
manner, these derivatives are also obtained at the streamline base 
point, point (5). Then, using the governing partial differential 
equations, the x-partial derivatives of the flow properties are found 
at points (1) to (5). For any bi characteristic which intersects the 
shock wave or the solid boundary, the time parameter [t (6 ) - t(k)] is 
found using equation (£.13) applied for i=l (i.e., the x-coordinate 
direction) while employing the appropriate intersection coordinates. 

At this stage, the system of compatibility equations may be solved 
for the flow properties at point (6) in a manner identical to that 
employed in the basic interior point scheme. 

The situaticn illustrated in Figure E.6 is quite general. In some 
instances, there are no bicharacteristic intersections with the solid 
boundary. Alternatively, there may be no intersections of the bi- 
characteri sties with the internal shock wave. There may be two bi char- 
acteristics intersecting with the shock wave, etc. 
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Another situation in which the interior point unit process must 
be applied in a modified form is illustrated in Figure E,7. In this 
figure, the Mach cone, with apex at the solution point. Intersects 
both the initial ‘Value plane and the internal shock wave. The point 
notation used in Figure E.7 is the same as that used in Figure E.6. 
However, in this case, the streamline base point, point (5), does not 
lie on the initial -value plane, but rather lies on the surface of 
the internal shock wave. 

The location of the streamline base point is obtained by extending 
the streamline from the initial-value plane to the surface of the shock 
wave. The point of intersection of the streamline with the shock wave 
is determined by employing the iterative scheme which is presented in 
Appendix D for finding a stream! ine-surface intersection point. That 
procedure requires that equation (E.68) be applied to three known 
shock wave solution points in the ini tial -value plane and three shock 
wave solution points in the current solution plane. Furthermore, 
interpolated values of the velocity components are required on the up- 
stream side of the shock wave at the point where the streamline inter- 
sects the shock wave. For this purpose, the following linear tri- 
variate interpolation polynomial is employed. 

f(x,y,z) = a 1 + a ? x + a^y + a 4 z (E.69) 

The coefficients a.. (i = l to 4) in equation (E.69) are determined 
by fitting this expression to four data points on the upstream side 
of the shock wave, as discussed in Appendix C. 

After the stream! ine-shock wave intersection point has been de- 
termined, the following fraction is formed 
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INITIAL- VALUE PLANE 



FIGURE E.7. SHOCK -MODIFIED INTERIOR POINT COMPUTATIONAL 

NETWORK (STREAMLINE BASE POINT ON SHOCK WAVE) 



(E. 70) 


e = [x s - x(5)]/(x s - Xj) 

where Xj and are tbe axial positions of the initial-value plane 
and the solution plane, respectively. If e is greater than a specified 
minimum value, an interior point unit process is performed on the 
downstream side of the shock wave. This unit process is almost iden» 
tical to that used for determining the solution at point (6) in 
Figure E.6. In this case, however, the streamline formula given by 
equation (E.12) is applied between the streamline-shock wave inter- 
section point and the solution plane. Interpolated flow property 
values at point {§) are determined by applying equation (E.67) to four- 
teen data points on the downstream side of the shock wave. 

If, on the other hand, e is less than the specified minimum value, 
an interior point unit process on the downstream side of the shock wave 
is not performed. Instead, the streamline from point (5) is projected 
onto the solution plane, and the flow properties at the solution point 
are determined by interpolation in the solution plane. The streamline 
integration from point (5) to point (6) employs equation (E.12). The 
flow property values at point (5) are obtained from equation (E.67) 
applied to fourteen, data points on the downstream side of the shock 
wave. Flow property values at the streamline-solution plane inter- 
section point are determined from the linear bivariate polynomial 

f(y,z) = a-j + a^y + a 3 z (E.71) 

The coefficients (i-1,2,3) in equation (E.71) are determined by 
fitting this expression to three data points in the current solution 
plane, as described in Appendix C. The order of integration for 
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determining the internal flow field is specified so that the downstream 
shock wave points and outer interior points in the downstream flow 
field sector are determined first. Tht; location of the solution point, 
in this case, is determined by an iterative loop which is terminated 
when the y and z coordinates of the projected solution point have 
converged . 


9. SHOCK-MODIFIED SOLID BOUNDARY POINT UNIT PROCESSES 

In some situations, the solid boundary point unit process must 
be applied in a modified form. One such application is illustrated in 
Figure E.8. In this situation, a portion of the Mach cone, with apex 
at the solid body solution point, intersects both the initial-value 
plane and the internal shock wave. The point notation used in Figure 
E.8 is identical to that used in Figure E.2, which depicts the computa- 
tional network for the standard body point unit process. The unit 
process employed in the present case is almost identical to the standard 
body point unit process. In the present case, however, the bi char- 
acteristic-shock wave intersection is handled in a manner identical to 
that employed in the shock-modified interior point unit process pre- 
sented in the previous section. 

In some situations, the entire Mach cone intersects the shock wave, 
as illustrated in Figure E.9. This situation occurs at a body point 
on the solution plane that is immediately downstream of a solid body- 
shock wave reflection, or at a body point on the solution plane that is 
immediately behind the shock wave emanating from the cowl lip. In the 
former case, the shock wave- sol id body intersection is a space curve 
in three-dimensions, whereas, in the latter case, the shock wave-solid 
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FIGURE E.8. SHOCK* MODIFIED SOLID BOUNDARY POINT 

COMPUTATIONAL NETWORK (TYPICAL APPLICATION) 




FIGURE E.9. SHOCK- MODIFIED SOLID BOUNDARY POINT 

COMPUTATIONAL NETWORK (POST SHOCK 
REFLECTION APPLICATION) 
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body intersection is a curve in a plane of constant x. The appropriate 
intersection algorithm is used as presented in Appendix D, and for the 
most part, procedures identical to those employed in the shock- 
modified interior point unit process are employed in this case. 

10. INTERNAL FLOW FIELD-SHOCK WAVE POINT UNIT PROCESSES 

Figure E.10 illustrates the overall computational network used in 
determining the solution for a typical shock wave point in the internal 
flow field. To determine the solution at the shock wave point, an 
interior point unit process must be performed to obtain the upstream 
flow properties at the location of the shock wave solution point. 

Figure E.1G illustrates both the computational network for the interior 
point unit process (denoted by primed numbers), and the computational 
network for the standard shock wave point unit process (denoted by 
unprimed numbers). The point notations employed in these computational 
networks are identical to those used in the corresponding standard 
unit processes. 

The computational procedure employed for determining the solution 
for an internal flow field-shock wave point is almost identical to the 
bow shock wave point unit process. The major difference between the 
two procedures is that for an internal flow shock wave point, the up- 
stream flow properties at the solution point are obtained from an 
interior point computation, rather than using free-stream data as in 
the bow shock wave point unit process. The required interior point 
unit process is essentially the same as the basic Interior point unit 
process presented in Section 4 of this appendix. In the present case, 
however, the streamline is not extended from a field point in the 
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FIGURE E.IO. INTERNAL FLOW FIELD-SHOCK WAVE POINT 

COMPUTATIONAL NETWORK (TYPICAL APPLICATION) 






initial -value plane to the solution plane, but rather it is extended 
from the shock wave solution point back to the initial-value plane. 

The position of the shock wave solution point is determined by the 
shock wave point unit process. To initiate the interior point computa- 
tion in the present case, flow property values are used from an 
adjacent field point in the flow field sector that is upstream of the 
shock wave in the solution plane. This modified interior point unit 
process requires searching the flow field sector upstream of the shock 
wave in the initial-value plane for the field point that is closest 
to the streamline-initial-value plane intersection point. This point 
is then used as the base point for the stencil of initial-value plane 
field points that are used in formulating the bivariate interpolation 
polynomial given by equation (E.11) (see Appendix C). 

For the first solution plane inside the inlet, the downstream 
bicharacteristic base point, point (1) in Figure E.ll, does not lie 
on the initial-value plane, but rather is located on the stream sur- 
face formed by the cowl boundary. To compute the pressure at point (2) 
from the wave surface compatibility relation, equation (E.59), the 
flow property values must be available at point (1), which requires 
that the flow property field must be known on the cowl surface. The 
body points on the cowl surface at the first internal flow solution 
plane, however, must be obtained from the unit process described in 
Section 9 of this appendix. That unit process requires that the 
flow property field on the downstream side of the shock wave be known. 
Hence, a simultaneous solid body point-shock wave point algorithm must 
be employed. This procedure was not developed in the present 
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FIGURE E.ll INTERNAL FLOW FIELD- SHOCK WAVE POINT 

COMPUTATIONAL NETWORK { FIRST INTERNAL FLOW 
SOLUTION PLANE APPLICATION) 


investigation. Rather, the shock wave points on the first internal 
flow solution plane are computed using a value of 41 in equation (E.37) 
equal to the value of 0 at the shock wave point in the initial- value 
plane which lies in the same meridional plane as the solution point. 

This provides a solution at each shock wave point on the first solu- 
tion plane without employing the compatibi li ty relation along the 
bicharacteristic. The body points on the cowl are then computed in the 
manner outlined in Section 9. On ensuing solution planes, except for 
the one immediately after a solid body-shock wave intersection, the 
bicharacteristic base point is located and the angle $ is iterated 
to convergence. 

When the internal shock wave intersects a solid boundary, as illus- 
trated in Figure E.12, a modification is required to the shock wave 
point unit process. In this case, instead of performing an interior 
point unit process to obtain the upstream flow properties at the solution 
point, a modified solid boundary point unit process must be employed. 
Moreover, the shock wave solution point, in this case, does not lie 
on the solution plane, but rather its position must be obtained by 
computing the intersection of the incident shock wave with the solid 
boundary. 

Finally, it should be noted that in order to achieve strict 
second-order accuracy in the internal flow shock wave point solution, 
global correction must be performed [this involves evaluating the 
cross derivatives at the solution point and using updated values of a. 
in equation (E.38)]. Time constraints in the present investigation 
did not permit the development of the global correction capability for 
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FIGURE E.I2. INTERNAL FLOW FIELD' SHOCK WAVE POINT 

COMPUTATIONAL NETWORK 

(INCIDENT WAVE 'BOUNDARY INTERSECTION APPLICATION) 
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the internal flow shock wave points. Hence, only local iteration can 
be performed for those points. 
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APPENDIX F 


OVERALL NUMERICAL ALGORITHM 


1 . INTRODUCTION 

The overall numerical algorithm consists of the repetitive 
application of the various unit processes to generate the global solu- 
tion for given boundary conditions and a specified set of initial 
data. 

The boundary conditions are represented by the fjrmulations pre- 
sented in Appendix D. The initial data are specified on a space-like 
plane of constant x. The x-coordinate axis is the longitudinal axis of 
the centerbody and the cowl. Moreover, the mean flow direction is 
assumed to be in the x-coordinate direction. 

An inverse marching scheme is employed in the overall numerical 
algorithm. The solution is obtained on space-like planes of constant 
x. The solution points on each plane represent the intersection 
points of continuous streamlines which are propagated from the data 
points specified on the initial ^iue plane. In addition to the 
streamline solution points, are the solution points representing the 
intersection of either the external or the internal shock wave with 
the solution plane. For the internal flow, the solution is also ob- 
tained on the space curves which represent the intersection of the 
internal shock wave with the solid boundaries. These space curves 
are defined by the locus of shock wave solution points. 
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Except in the vicinity of a shock wave reflection with a solid 
boundary, the axial (x) distance between successive solution planes 
is determined by the application of the Courant-Friedrichs-Lewy 
(CFL) stability criterion. In the vicinity of a shock wave intersec- 
tion with a solid boundary, the axial step is controlled by special 
constraints which insure that the entire shock wave-solid boundary 
intersection falls between two adjacent solution planes. 

After each solution plane is computed, the mass flow rate across 
that plane is calculated using trapezoidal rule integration. Constancy 
of the overall mass flow rate in the internal flow field computation 
gives an indication of the overall accuracy of the numerical integra- 
tion. The stagnation pressure and stagnation temperature are 
calculated at each solution point. For the adiabatic flow of a 
calorically perfect gas, the stagnation temperature should remain 
constant. 

In the numerical analysis, the flow field is divided into two 
regimes: the internal flow regime and the external flow regime, as 

illustrated in Figure F.l. The flow field integration in each of these 
two regimes is controlled by separate logic modules in the computer 
program. Toe forebody flow field integration is performed first. 

Then, the internal flow field is computed. The computer program 
developed in the present investigation has the capability to perform 
the internal flow field integration with or without the discrete 
fitting of the internal shock wave system. The option in which shock 
waves are not discretely fitted might be employed if the internal 
shock waves are of relatively weak strength, and thereby an acceptable 
solution could be obtained by smearing the internal discontinuities. 





FIGURE F.l. MIXED- COMPRESSION AIRCRAFT INLET 


From a computation point of view, the internal flow field in 
which shock waves are not discretely fitted is the easiest solution 
to compute. For flow fields in which shock waves are discretely 
fitted, the external flow about the forebody is less difficult to ob- 
tain than the internal flow, since in the external flow the shock wave 
represents a bound to the computational regime. Discrete fitting 
of the shock wave throughout the computational regime, as is done in 
the internal flow field integration, greatly complicates the numerical 
algorithm. 

In this appendix, the overall control logic used in each of the 
three flow field integration options is discussed. Regulation of the 
axial marching step size, generation of the initial-value surface 
data, and considerations of flow symmetry are also discussed. All of 
the unit processes referred to in this appendix are discussed in 
Appendix E. 

2. C0URANT-FRIEDR1CHS-LEWY (CFL ) STABILITY CRITERION 

Except in the vicinity of an internal shock wave-solid boundary 
intersection, the axial marching step between successive solution 
planes is determined by the application of the Courant-Friedrichs-Lewy 
(CFL) stability criterion (9), The CFL stability criterion will be 
satisfied at each solution point if the convex hull of the finite 
difference network contains the differential zone of dependence of 
the solution point. The convex hull of the finite difference network, 
illustrated in Figure F.2, is defined by the outer periphery of 
initial-value plane field points used in determining the fit point 
stencil for the quadratic bivariate interpolation polynomial. The 
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FIGURE F.2. CFL STABILITY CRITERION 


differential zone of dependence, also illustrated in Figure F.2, is 
the region defined by the intersection of the Mach :one (whose apex 
is at the solution point) with the initial-value plane. 

The maximum allowable marching step for each streamline is the 
x-step for which the Mach cone just touches the convex hull. That 
step size is given by 

Ax = [u 2 /(cq)][l - (c/q)(q 2 /u 2 - l) 1/2 ]R min (*M) 

where Ax is the maximum allowable axial step, u is the x-component 
of the velocity, q is the velocity magnitude, and c is given by 

c 2 = a 2 q 2 /(q 2 - a 2 ) (F.Z) 

where a is the local sonic speed. In equation (F.l), R^ n is the 
distance from the streamline base point in the initial-value plane 
to the nearest field point on the convex hull of the finite differ- 
ence network (see figure F.2). 

Equation (F.l) is applied at every streamline solution point, 
the actual marching step being selected as the ax value at the most 
restrictive point. It should be noted that this expression is applied 
only to streamline points, the shock wave points being excluded. 
Furthermore, in the internal flow field integration, the shock wave 
points are ignored in defining the convex hull of the finite differ- 
ence network when application of the stability criterion is made to 
a streamline point. 
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3. INITIAL- VALUE PLANE 


The initial data are specified on a plane of constant x. The 
flow must be supersonic at every point on this plane. For uniqueness 
and existence of a genuine solution* the values of the five dependent 
variables (u» v, w, P, and p) prescribed on this surface must have 
at least continuous first derivatives. 

If the forebody flow field is to be determined, the initial -value 
plane must be specified at an axial (x) station that is upstream of 
the forebody computational flow regime (see Figure F.l). The 

solution is then found along the streamlines that pass through 
the data points specified on the initial-value plane, although some 
streamline addition and deletion are performed on the ensuing solu- 
tion planes as described in Section 5 of this appendix. 

If only the internal flow field is to be determined, the initial- 
value plane must be specified at the axial station which corresponds 
to the x- position of the cowl lip (see Figure F.l). The cowl lip is 
assumed to be contained in a plane of constant x. For the integration 
of the internal flow field, a point redistribution is performed on 
the initial-value plane. This point redistribution is required in 
order to have streamlines which lie in the stream surface formed by 
the cowl boundary. The solution is then found along the streamlines 
that pass through the redistributed points on the plane at the cowl 
lip axial station. 

The initial-value plane may be specified by the user, or if the 
forebody is conical up to the axial station where the initial -value 
plane is located, the flow property field on the initial -value plane 
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can be generated internally in the computer program. The internally 
generated initial -'.slue plane is obtained by an approximate technique 
which employs the Taylor-Maccoll solution for the flow about a circu- 
lar cone at zero incidence, A superposition procedure is used to ob- 
tain an approximation to the flow about a circular cone at nonzero 
angle of attack by neglecting the cross flow effects. This superposi- 
tion procedure effectively amounts to computing the flow turning 
angle in the meridional plane of the given solution point, and then 
obtaining the flow properties at that point by applying the Taylor- 
Maccoll solution for a cone half-angle equal to the flow turning 
angle. The shock wave angle is then measured from the original 
streamline direction in the appropriate meridional plane. It must be 
emphasized that this is only an approximate technique, giving the well 
accepted Taylor-Maccoll solution at zero incidence, but becoming 
increasingly less accurate as the angle of attack is increased. 

The solution obtained by Jones (33) for the flow about a circu- 
lar cone at nonzero incidence has been well substantiated. Using a 
conversion algorithm, the results of the computer program developed by 
Jones can be made compatible with the input data required by the 
computer program developed in the present investigation. Many of the 
computed results presented in Section VI were obtained using the 
results of Jones' program as initial data. For situations ir> which 
the forebody is conical up to the axial station where the initial- 
value plane *is located, the Jones program is the recommended source 
for the initial data. 
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If the forebody Is not conical ahead of the axial station of the 
initial-value plane, another source of initial data must be used. 

If available, experimental data may be employed. 

4. SOLUTION PLANE POINT NETWORK AND FLOW SYMMETRY 

The computational point network is based on a series of circum- 
ferential and radial stations. The point networks for the various 
flow symmetry options are illustrated in Figure F.3. In this figure, 
the index i corresponds to the ith circumferential station and the 
index j corresponds to the jth radial station. In all cases, the 
streamlines on the surface of the center body are denoted by j = 1. 


v 


For the forebody flow field, the bow shock wave solution points are 
denoted by j = n. For the internal flow field, the streamlines on 
the surface of the cowl are denoted by j = n. The computed sector, 
in general, is bounded by the circumferential stations corresponding 
to i = 1 and i = m. This point arrangement produces a rectangular 
logic array in the computer program. 

The points at any circumferential station in axisymmetric flow, 
or on a plane of flow symmetry in three-dimensional flow, lie on a 
straight line. Moreover, for axisymmetric flow, the radial stations 
correspond to circular rings. In general, however, the solution 
points at a given circumferential station do not lie on a ray, nor 
do the radial stations correspond to circular rings. 

For the internal flow option in which shock waves are discretely 
fitted, the shock wave solution points are also represented in this 
point arrangement. Special logic is used in the computer program 
such that the shock wave solution points float in the storage arrays 


209 





210 


(X-AXIS OUT OF PAGE) 



(a) NO PLANES OF SYMMETRY (b) ONE PLANE OF SYMMETRY 


FIGURE F.3. COMPUTATIONAL POINT NETWORKS 














as the shock wave travels between the centerbody and cowl on successive 
solution planes. Gn a given solution plane, the shock, wave solution 
points at adjacent circumferential stations do not, in general, have 
to lie at the same radial station. 

The computer program takes advantage of flow symmetry when it 
exists in the flow field. In these instances, the entire solution 
plane does not have to be computed, but rather only an appropriate 
section of it. The remaining sections of the solution plane may be 
obtained by reflection of the points in the computed sector. This 
procedure yields a significant reduction in computer execution time. 

The four flow symmetry options that have been incorporated into 
the analysis are depicted in Figure F.3. Figure F.3(a) illustrates 
the most general case when no flow syrtmetry is present. Figure F . 3 ( b ) 
illustrates the case when one plane of flow symmetry is present. In 
this case the computed sector is the half-plane bounded by the 
y-axis and containing the +z-axis. The integration region in this 
case is bounded by the i = 1 circumferential station on the 
+y-axis and by the i = m circumferential station on the -y-axis. 

This case of flow symmetry is the one most likely to arise in the 
class of problems being considered in this investigation. Figure 
F.3(c) illustrates the case when two planes of flow symmetry are pre- 
sent. This option would be used to compute the flow field about asym- 
metric bodies at zero angle of attack. In this instance, the computed 
sector is the quadrant bounded by the +y-axis and the +z-axis. The 
circumferential station corresponding to i = 1 lies on the +y-axis 
and the circumferential station corresponding to i = m lies on the 
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+z-axis. Finally, Figure F. 3(d) illustrates the axisynmetrie flow 
option where the computed sector is limited to the single circumfer- 
ential station (ray) lying on the +y-axis. This option would be used 
to compute the flow field about axi symmetric bodies at zero angle of 
attack. 

The numerical algorithm does not apply special unit processes 
when a solution point lies on a plane of symmetry. Rather, a point 
reflection about the plane of symmetry is performed in the initial- 
value plane, and the appropriate unit process is then applied in 
standard form. This procedure yields satisfactory results and elimi- 
nates the need for devising special unit processes. 

5. EXTERNAL FLOW ABOUT THE FOREBODY 

With the forebody geometry specified and the flow property field 
on the initial-value plane determined, the external flow about the 
forebody can be calculated. In the computation of this flow field, 
the distance between successive solution planes is determined by the 
application of the CFL stability criterion. The last solution plane 
in the forebody flow field computation is made to coincide with the 
x-position of the cowl lip. 

After the axial step between the current initial-value plane and 
the current solution plane has been determined, the solid boundary 
point unit process (see Appendix E) and the interior point unit 
process (see Appendix E) are applied. These unit processes achieve 
second-order accuracy without the need for global iteration. Hence, 
these unit processes are applied at the appropriate points until con- 
vergence is obtained without using information from neighboring points 
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in the solution plane. 

Once the solution at each solid boundary point and interior point 
has been determined, the bow shock wave point unit process (see 
Appendix E) is applied at each shock wave solution point in the com- 
puted sector. Global correction is then applied for these points, 
if desired. The position of each shock wave solution point is made 
to lie in the meridional plane defined by the outer-most interior 
field point which is on the same circumferential station as the shock 
wave point. As a consequence, in axi symmetric flow, the streamline 
and shock wave solution points on a given circumferential station lie 
in the same meridional plane on all succeeding solution planes. In 
three-dimensional flow, however, except on a plane of flow symmetry, 
the solution points corresponding to a given circumferential station 
do not lie in the same meridional plane on successive solution p^nes. 

!n the forebody flow field integration, periodic streamline 
addition and deletion are performed. The streamline addition is 
required to retain a well -dispersed computational mesh, since at suc- 
cessive solution planes more and more mass is captured. Moreover, 
convergence of the streamlines towards the forebody occurs as the 
flow progresses downstream. Periodic point deletion is required since 
the continued addition of streamlines would produce an excessively 
large number of computational mesh points, thereby unduly increasing 
computer execution time and machine storage requirements. The stream- 
line addition and deletion procedures are outlined in the following. 

A depiction of a typical forebody flow streamline pattern is given in 
Figure F.4. 
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For the purposes of point addition, after the points on the solu- 
tion plane have been computed, the mass flow rate across that plane 
is calculated. If this mass flow rate is significantly larger than 
the mass flow rate across the last solution plane where point redis- 
tribution was performed, a new ring of solution points is added 
between the ring of shock wave solution points (j = n) and the ring 
of outermost interior field solution points (j = n - 1). The co- 
ordinates of each of these inserted solution points is obtained by 
forming the arithmetic average of the coordinates of the shock wave 
solution point and the outermost interior field point corresponding 
to the circumferential station of the new point. The flow properties 
at each of the inserted solution points are obtained by interpolation 
using the quadratic bivariate polynomial 

f(y.z) = a 1 + a^y + a^z + a 4 yz + a § y 2 + a g z 2 (F.3) 

where f(y,z) denotes a general function of the coordinates y and z. 

The coefficients a. (i=l to 6) in equation (F.3) are obtained by a 
least squares fit of nine data points in the solution plane, as de- 
scribed in Appendix C. 

Point deletion occurs when the number of radial stations has 
reached a specified limit. In point deletion, the body streamline 
points are retained in storage, while selected interior streamline 
points are deleted from storage. Refinement of this technique is pro- 
vided by having two limits to the number of allowable radial stations. 
The first limit is employed when the mass flow rate at the given 
solution plane is less than a specified fraction of the estimated 
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flow rate at the cowl lip. The second and larger limit is employed 
when that fraction has been exceeded. 

Finally, it should be noted that the influence of molecular 
diffusion can be included in the forebody flow field computation. 

6. INTERNAL FLOW IN WHICH SHOCK WAVES ARE NOT DISCRETELY FITTED 

The program option in which the internal flow field is computed 
without the discrete fitting of the internal shock wave system might 
be employed in the cases where the internal shock waves are weak in 
strength, and thereby an acceptable solution could be obtained by 
smearing all internal discontinuities. This option requires that 
only two unit processes be employed: the interior point unit process 
and the solid boundary point unit process. The influence of molecular 
transport can be included in the computation of this flow field. 

The initial -value plane of the internal flow computation is 
constituted by the last solution plane of the forebody flow field 
integration. Alternatively, the initial-value plane may be specified 
at the cowl lip axial station without employing the forebody flow 
field integration option. This technique is recommended if the fore- 
body is conical up to the cowl lip axial station. 

The computer program developed in the present investigation assumes 
that the bow shock wave falls outside of the cowl lip, or, in the 
limit, intersects the cowl lip exactly. The program does not have the 
capability to compute the internal flow field when the bow shock wave 
has been ingested into the annulus. 

With the initial -value plane specified, a point redistribution 
on this plane is performed to obtain a uniform point distribution and 

217 



i 

i 

i 

i 

| 

1 




f 


jusst «* 


-flftfGttrt; ■a««W6i *aiW5fcri ‘ *»*«*■*• --tniftm '.■****& jwjku tsevc'i.-i *' » - sKiae** « ,v-WiM' -?** v - * ^4 Jgwi “■*****•<* ..'^jwt .^rar**_ , 


to obtain streamlines which lie in the stream surface formed by the 
cowl boundary. The redistributed points are arranged symmetrically 
in the computed sector. These points lie on rays which have equal 
angular increments from one another, with the points on each ray 
being spaced at equal radial increments. The radial station j = 1 
corresponds to the centerbody streamline points, and j - n corresponds 
to the cowl streamline points. The properties at these points are 
obtained by interpolation. 

With the point redistribution performed, the internal flow field 
integration proceeds in a manner similar to the external flow field 
integration, except that only two unit processes are used: the 

interior point unit process and the solid boundary point unit process. 
No point addition or deletion is performed. The internal flow field 
integration is terminated either when a specified axial station is 
reached or when the flow becomes subsonic, 

7. INTERNAL FLOW IN WHICH SHOCK WAVES ARE DISCRETELY FITTED 

A point redistribution is first performed on the initial-value 
plane at the axial station of the cowl lip as described in the previous 
section. After the upstream flow properties have been determined at 
each of the cowl lip solution points in the computed sector, the 
downstream flow properties are obtained at each of these points by 
application of the solid body-shock wave point unit process. 

In the integration of the internal flow field in which shock waves 
are discretely fitted, the axial step is obtained by the application 
of the CFL stability criterion, except in the vicinity of a shock 
wave reflection, where special constraints are employed. After the 
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axial station of the solution plane has been determined, the internal 
shock wave is projected from the current initial -value plane to the 
current solution plane in the meridional plane passing through the 
x-axis and the previous shock wave point on the initial -value plane 
as illustrated in Figure F.5. The location of the shock wave solu- 
tion point is obtained by applying the following equation. 

dr s /dx = tanSj (F.4) 

In equation (F.4), dr^ is the increment in radius between the pro- 
jected shock wave point and the previous shock wave point on the 
initial-value plane, dx is the corresponding increment in axial dis- 
tance, and gj is the angle subtended by the shock wave and the x-axis 
at the initial -value plane shock wave point and in the meridional 
plane defined by the ini tial -value plane shock wave point. Equation 
(F.4) is applied for each shock wave point in the computed sector, 
thereby yielding the locus of projected shock wave points in the 
solution plane. Interpolated values of the shock wave radius in the 
solution plane are obtained by employing the following equation. 

r s (0) = a 1 + a 2 e + a 3 © 2 (F.5) 

In equation (F.5), r s (0) is the shock wave radius at the polar angle 
0 = tan"^(z/y), and the coefficients a- (i=l,2,3) are obtained by 
fitting this expression to three projected shock wave points, as de- 
scribed in Appendix C. Equation (F.5) is applied at every circumfer- 
ential station in the computed sector. Hence, the shock wave location 
in the solution plane is represented by a series of overlapping one- 
dimensional curve fits. 
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After the tentative position of the shock wave in the solution 
plane has been determined, the streamlines that are in the flow field 
sector that is upstream of the shock wave in the initial -value plane 
are projected from the initial-value plane to the solution plane, as 
illustrated by streamlines 1 to 6 in Figure F.5(a) and by streamlines 
9 to 13 in Figure F. 5(b) . This is accomplished by applying the equa- 
tion of a streamline 

dx. = u.dt { l ~ 1 ,2,3) (F.6) 

where (1=1 ,2,3) denotes the three cartesian coordinates x, y, and 
z, respectively, u^. (i=l,2,3) denotes the corresponding velocity com- 
ponents u, v, and w, respectively, and t is the time of travel of a 
fluid particle along the streamline. Equation (F.6) is first applied 
in the x direction. Since the axial step dx is known from the appli- 
cation of the CFL stability criterion, the time parameter dt may be 
determined. Then, application of equation (F.6) for the y and z 

directions allows the y and z coordinates of the projected streamline 

2 2 1/2 

point to be computed. The radius r = (y + z ) ' and polar angle 
9 = tan’^fz/y) of each of the projected streamline points are then 
computed. 

The radius of the projected streamline point is then compared to 
the radius of the shock wave, given by equation (F.5), in the meri- 
dional plane defined by the projected streamline point. If the pro- 
jected streamline point is in the upstream flow field sector on the 
solution plane (i.e., the streamline does not intersect *he shock 
wave), then a standard interior point or solid boundary point unit 
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process is applied to obtain the solution at this point. Jf the 
streamline appears to intersect the shock wave, as illustrated by 
streamlines 5 and 6 in Figure F.5(a) and streamlines 9 and 1G in 
Figure F.5(b), then the application of the unit process to determine 
the solution is deferred. 

At this stage, the upstream and downstream shock wave solution 
points are determined at each circumferential station in the solution 
plane computed sector using the internal shock wave point unit process. 
This procedure defines the property field on both the upstream and 
downstream sides of the internal shock wave. 

Next, the body streamline solution points are computed at every 
circumferential station in the downstream flow field sector on the 
solution plane. In some instances, computing the solution at these 
points may entail using flow property information from the downstream 
side of the internal shock wave if the Mach cone, with apex at the 
solution point, intersects the shock wave surface. Determining the 
solution at each of these points thereby defines the flow property 
field on the boundary stream surface in the downstream flow field 
sector. 

At l..is stage, the solution on each of the streamlines which have 
not yet been computed is determined. The streamlines that are in the 
downstream flow field sector on the initial-value plane will remain 
in the downstream flow field sector on the solution plane (see Figure 
F.5). The solution at these points is determined by the application 
of the standard interior point unit process, unless a portion of the 
Mach cone, with apex at the solution point, intersects the internal 
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shock wave or the solid boundary, in which case the modified interior 
point unit process is applied. For streamlines which penetrate the 
internal shock wave [streamlines 5 and 6 in Figure F. 5{a ) and stream- 
lines 9 and 10 in Figure F.5(b)], the appropriate modified interior 
point unit process is applied. For the streamlines whose solution 
was deferred due to a possible shock wave penetration, but which 
ultimately did not intersect the shock wave, the standard interior 
point scheme is applied. The solution points are ordered in the 
storage arrays in the order of increasing radius on a given circum- 
ferential station. So a post computation interchange of the stream- 
line solution points with the shock wave solution points is performed 
for the streamlines which initially appeared to intersect the shock 
wave but ultimately did not. 

The process just outlined is applied repetitively until the 
internal shock wave intersects a solid boundary. Special logic is 
used in the computer program for the computation of a shock wave 
reflection. The overall scheme used in this case is now presented. 

The initial step in the computation of the shock wave-solid 
boundary reflection is to obtain an estimate of the axial location, 
at a discrete number of points, where the incident shock wave inter- 
sects the solid boundary. Except for the case of ax i symmetric flow, 
the intersection of the incident shock wave with the solid boundary 
defines a three-dimensional space curve, as illustrated in Figure F.6. 
In axisymmetric flow, this curve lies in a plane of constant x. 

Points along the space curve are determined by obtaining the inter- 
section of the shock wave and the solid boundary, where both of these 
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surfaces are represented as straight line segments in the meridional 
planes passing through the shock wave points in the initial -value 
plane. For a given meridional plane, the shock wave is represented 
by equation (F.4), where dr £ is the increment in radius between the 
shock wave-body intersection point and the shock wave point in the 
initial-value plane, dx is the corresponding increment in axial 
distance, and gj is the angle subtended by the shock wave and the x- 
axis in the meridional plane defined by the appropriate shock wave 
solution point in the initial -value plane. The local body surface is 
approximated in the meridional plane by the equation 

dr^/dx * m (F . 7 ) 

where dr^ is the change in the radius of the body between the shock 
wave-body intersection point and the body point in the initial-value 
plane, dx is the corresponding increment in axial distance, and m is 
the local slope of the body in the given meridional plane. Equations 
(F.4) and (F.7) are solved simultaneously to obtain the intersection 
point in the given meridional plane. The intersection point for 
every meridional plane defined by the shock wave points on the initial- 
value plane is so determined. The locus of these intersection points 
determines the space curve illustrated in Figure F.6. 

At this stage, the points on the space curve which are nearest 
to and farthest away from the initial-value plane are determined. If 
the axial distance between the nearest point and the initial-value 
plane is greater than a specified fraction of the marching step 
allowed by the CFL stability criterion, then another solution plane 


is computed, the location of this plane being just slightly upstream 
of the shock wave-body intersection. The entire procedure outlined 
above is then repeated. Alternatively, if the distance between the 
nearest shock wave-body intersection point and the initial -value plane 
is less than this fraction of the allowable marching step, then the 
axial position of the next solution plane is selected such that the 
space curve representing the incident shock wave- body intersection 
is entirely contained between the initial-value plane and the solution 
plane. At high angles of attack, this procedure may require that the 
axial step between the initial-value plane and the solution plane 
be greater than that allowed by the CFL stability criterion. This 
implies that the Courant number, which is the ratio of the axial step 
taken to the axial step allowed by the CFL stability criterion, is 
greater than unity. To maintain an effective Courant number less 
than unity, the fit point stencils used in the univariate, bivariate, 
and trivariate interpolation polynomials are adjusted in accord with 
the Courant number of the actual step taken. That is, if the Courant 
number is approximately two, then every other point is used in the 
interpolation fit point stencils instead of the inmediate neighbors 
(which correspond to a unity Courant number), etc. This ensures that 
the convex hull of the finite difference network engulfs the differ- 
ential domain of dependence, thereby satisfying the CFL stability 
criterion. 

After the axial position of the solution plane has been determined 
and the Courant number computed, the internal shock wave point unit 
process is applied at every circumferential station in the computed 


227 


i««aaar.tfgpwrn ■Tttmw.sEBBartaaaBaw "sseoa jawi — 


sector at the Intersection of the incident shock wave with the solid 
boundary. This procedure defines the property field on both the up- 
stream and downstream sides of the incident shock wave. 

At this stage, the initial-value plane upstream sector body 
streamlines are extended from the initial -value plane to the space 
curve defined by the intersection of the incident shock wave with the 
solid boundary, as illustrated in figure F.6. The solution for both 
the upstream and downstream shock wave properties has been obtained 
on the space curve by the application of the internal shock wave 
point unit process. Hence, both the upstream and downstream properties 
at the points where the body streamlines intersect the space curve 
may be found by interpolation, for this purpose the following quad- 
ratic univariate polynomial is employed 

f (6 ) = a ] + a ? o + a 3 e ? (F.8) 

where f {© ) denotes a general function of the polar angle 8. The 
coefficients a.. (i-1,2,3) in equation (F.8) are obtained by fitting 
this expression to three data points on the space curve as described in 
Appendix C. To determine the intersection point of the body streamline 
with the space curve, an iterative technique is used. Moreover, after 
each iteration, the projected streamline point is adjusted along the 
direction of the body normal projection in the (y,z)-plane such that 
the streamline point lies on the boundary surface. Equation (F.8) 
is applied for both the upstream and downstream shock wave properties. 
Hence, the incident shock wave downstream properties are known at the 
body streamline points. 


228 



>sf.t ^simanmm s'sees m ■ fluastSSa 


At this stage, the solid body-shock wave point unit process is 
applied at each of the body streamline points in the computed sector 
that are on the space curve. This defines the reflected shock wave 
downstream properties at the body streamline points on the space 
curve . 

Using a procedure similar to that used previously, the shock wave 
is then projected from the space curve to the current solution plane. 
This projection is performed in the meridional planes containing the 
body streamline points on the space curve. This procedure yields the 
tentative shock wave shape in the solution plane- 


At this stage, the body streamline points in the solution plane 
that are in the downstream flow field sector in the initial-value 
plane are computed by use of the solid boundary point unit process 
(see Figure F.7). This unit process is applied at every such point 
in the computed sector. As a consequence, the flow property field on 
the stream surface formed by the solid boundary is defined. 

Next, the remaining streamlines that are in the initial-value 
plane downstream flow field sector are projected from the initial-value 
plane onto the solution plane. A test is then made to determine 
whether or not each of these streamlines intersects the reflected 
shock wave (see Figure F.7). Those streamlines which do not intersect 
the reflected shock wave will lie in the upstream flow field sector 
on the solution plane (points 7 to 13 in Figure F.7). The solution 
at these points is determined using the standard interior point 
scheme, or if the Mach cone, with apex at the solution point, inter- 
sects the incident shock wave or solid boundary, the appropriate 
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FIGURE F.7. TYPICAL STREAMLINE NETWORK AT AN 

INTERNAL SHOCK WAVE REFLECTION 







modified interior point unit process is applied. Those streamlines 
which appear to intersect the reflected shock wave have their compu- 
tation deferred. 

At this stage, the upstream and downstream shock wave points are 
determined at every circumferential station in the solution plane com- 
puted sector. This procedure defines the property field on both the 
upstream and downstream sides of the reflected internal shock wave. 

Next, the solution is obtained at each body streamline point in 
the downstream flow field sector on the solution plane {see Figure 
F.7). The modified solid boundary point unit process is applied in 
this situation, which requires using flow property information on the 
downstream side of the reflected shock wave. After the application 
of the body point unit process at each point in the computed sector, 
the property field on the solid boundary is defined. 

At this stage, the streamlines that are in the downstream flow 
field rector in the initial-value plane and that intersect the re- 
flected shock wave are computed. These points require using the 
modified interior point unit process and use flow property information 
on both the upstream and downstream sides of the reflected internal 
shock wave (see Figure F.7). 

Finally, the streamlines that are in the upstream flow field 
sector in the initial-value plane are extended to the surface of the 
incident shock wave and their respective intersection points with this 
surface are determined (see Figures F.7 and F.8). These streamlines 
are then extended from the downstream side of the incident shock wave 
to the current solution plane. If the projected streamline does not 
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intersect the reflected shock wave, a modified interior point unit 
process is applied using flow property information on the downstream 
side of the incident shock wave. If the projected streamline inter- 
sects the reflected shock wave, the intersection point is found with 
this surface. A modified interior point unit process is then applied 
on the downstream side of the reflected shock wave. 

After all of the points have been determined on the solution 
plane that is immediately downstrear of the shock wave-solid body 
reflection, control is returned to the driving algorithm until another 
shock wave-solid body reflection is encountered. 

Figures F.6 to F.8 illustrate the intersection of the shock wave 
with the eenterbody. Similar results hold when the shock wave inter- 
sects the cowl. 

The internal flow field integration is terminated when either a 
specified axial station is reached or when the flow becomes subsonic. 
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APPENDIX G 

CALCULATION OF THE TRANSPORT TERMS 


1 . INTRODUCTION 

The numerical procedure developed in this investigation has the 
capability to include the influence of molecular transport on the 
solution by treating the viscous and thermal diffusion terms in the 
governing equations as forcing functions, or correction terms, in the 
method of characteristics scheme. At present, the computer program 
has the capability to include the influence of viscous and thermal 
diffusion in the computation of the external flow field about the 
forebody, and in the computation of the internal flow field in which 
shock waves are not discretely traced. The program option which per- 
forms discrete fitting of the internal shock wave system does not have 
the capability to include the influence of molecular transport in the 
computation, but rather assumes the flow to be inviscid and adiabatic. 


2. EXPRESSIONS FOR THE TRANSPORT TERMS 

The expressions for the transport forcing functions are derived in 
Appendix A, and are summarized below. 


F * = 
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where 


r = 1_ (<£) 

^ pT ; 3s ; p 


(6.5) 


In equations (G.l) to (G.5), u, v, and w denote the velocity components 
in the x, y, and z coordinate directions, respectively, P is the pres- 
sure, p denotes the density, T is the absolute temperature, s denotes 
the entropy per unit mass, u represents the dynamic viscosity, and k is 
the thermal conductivity. The subscripts x, y, and z on the right-hand 
sides of equations (G.l) to (G.4) denote partial differentiation in the 

corresponding coordinate direction, whereas F , F , and F on the left- 

x y z 

hand sides denote the transport forcing functions in the x, y, and 
z component momentum equations, respectively. F g is the transport 
forcing function in the energy equation. 
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3* COMPUTATION OF THE TRANSPORT FORCING FUNCTIONS 

During the course of the program development* a number of methods 
were deviled in an effort to obtain a good approximation to the trans- 
port forcing functions. One such method was based on employing a 
quadratic trivariate interpolation polynomial whose coefficients were 
determined by a least squares fitting of a number of known field points 
on the initial -value plane and the previous solution planes. This 
polynomial was employed to determine the five dependent properties u, 
v, w, P, and p. The spatial derivatives of the velocity components ap- 
pearing in equations (G.l) to (6.4) were then obtained by analytically 
differentiating the respective interpolation polynomials. Spatial 
gradients of pressure and density were obtained in a similar manner. 
Then, by differentiation of the thermal equation state, temperature 
derivatives were expressed in terms of the pressure and density deriva- 
tives. The molecular transport properties and their spatial gradients 
were obtained using the procedures presented in Appendix A. 

This method of calculating the transport forcing terms was con- 
sidered to have good accuracy. The computer execution time required by 
this method, however, was felt to be unacceptable. This prohibitive 
execution time was primarily due to the least squares curve fitting of 
the tri variate interpolation polynomials. Consequently, a more 
efficient method with acceptable accuracy was sought for approximating 
the transport terms. The method which was selected is presented below. 

For the interior point and solid boundary point unit 
processes, the transport terms must be computed at all points in the 
computational network (see Figures E.l and E.2). For the bow shock 
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wave point unit process, the transport terms must be computed at the 
solution point and at the intersection point of the bi characteristic 
with the initial-value plane (see Figure E.3). For each of these 
unit processes, partial derivatives of the dependent properties with 
respect to y and z in the initial-value plane are obtained by analyti- 
cally differentiating the quadratic bivariate interpolation polynomial 

f(y»z} = ^ + a^ + a 3 z + a 4 yz + a g y 2 + a g z 2 (G.6) 

The coefficients a^ (1=1 to 6) in equation (G.6) are determined by a 
least squares fit of nine data points in the initial-value plane as 
discussed in Appendix G. Equation (G.6) is applied for the five 
dependent flow properties u, v, w, P, and p. Spatial derivatives of 
pressure and density are required [even though they do not appear 
explicitly in equations (G.l) to (G. 4)] because spatial derivatives 
of temperature are expressed in terms of pressure and density deriva- 
tives through differentiation of the thermal equation of state as 
discussed in Appendix A. 

In the solution plane, partial derivatives of the dependent 
properties with respect to y and z are equated to the corresponding 
derivatives in the initial-value plane. For the interior point and 
boundary point schemes, the derivatives at the solution point are set 
equal to those at the streamline base point. For the bow shock wave 
point scheme, the solution point derivatives are equated to those at 
the bicharacteristic base point. The evaluation of these derivatives 
in the solution plane would require that a global iteration algorithm 
be employed. In this algorithm, the property field on the solution 
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plane would first be determined by a predictor application of the 
appropriate unit process at each point in the computed sector. Then, 
by fitting equation (6.6) to solution plane field points the appropriate 
partial derivatives could be obtained. In a similar manner, ensuing 
corrector applications would be performed until overall convergence was 
achieved. The attendant increase in algorithm complexity and computer 
execution time using this global iteration procedure, however, was felt 
to be unwarranted since the transport terms are assumed to be of 
secondary importance in determining the solution. 

Partial derivatives with respect to x in equations (6.1) to 
(6.4) are obtained from the following quadratic univariate interpolation 
polynomial . 

f (x ) = a^ + a 2 x + a 3 x 2 (6.7) 

The coefficients a^ (1=1, 2, 3) in equation (6.7) are determined by fitting 
this expression to three data points. The first data point is located 
on the solution plane that is immediately upstream of the current 
initial-value plane, the second data point is on the initial-value 
plane, and the third data point is the solution point itself. For the 
interior point and boundary point unit processes, the fit points are 
located on the streamline which passes through the solution point. For 
the bow shock wave point unit process, the fit points are the shock 
wave solution points corresponding to the circumferential index of the 
solution point. Special logic in the computer program takes account 
of point deletion and addition in the forebody flow field computation 
and thereby insures that the appropriate fit points are selected. Of 
course, for a predictor application of either the interior point or 
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boundary point unit processes, property values at the solution point 
are equated to those at the streamline base point in the initial-value 
plane. 

Equation (G.7) is applied for the five dependent flow properties 
u» v, w, P, and p. Analytical differentiation of equation (G.7) yields 
approximations to the x-partial derivatives. Differentiation of the 
thermal equation of state allows the spatial derivatives of temperature 
in the x -coordinate direction to be expressed in terms of the correspond 
ing pressure and density derivatives. This formulation yields an x- 
partial derivative which is constant in a given x-plane. 

Since equation (G.7) uses data on a previous solution plane, 
derivatives cannot be evaluted using this representation until at least 
one previous solution plane is available. Furthermore, the derivatives 
obtained using this formulation are only approximations to the x-partial 
derivatives since the y and z coordinates of each of the three fit 
points are not, in general, identical. Considering that the effects of 
molecular diffusion are assumed to be small, this approximation is 
acceptable. 

The molecular transport properties and their spatial gradients 
are obtained using the procedures presented in Appendix A. 
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APPENDIX H 


NOMENCLATURE 


ENGLISH SYMBOLS 


3^ »b - tC . i d.. 


B, 


c 

e 


f 

F 


,F ,F 

x y 


z 


* 


F 


e 


i,j.k 


i' 



l 


t ~ ( A ^ » l i , £ ^ * ) 


M 


n = ( n x * n y » n z ) 
n i 

b = (n bx’ n by ,n bz* 


n 


bi 


sonic speed 

general curve fit coefficients 
body force vector in index notation 
velocity of divergence of Mach conoid surface 
internal energy per unit mass 
general interpolation polynomial function 
forcing functions in the x, y, and z component 
momentum equations and energy equation, 
respectively 

unit vectors in the x, y, and z directions, 
respectively 

unit vectors in the x', y ' , and z' directions, 
respectively 

unit vector along the space curve defined by 
the intersection of the shock wave with either 
the initial -value plane or a solid boundary 
Mach number 

unit vector normal to a wave surface 
above unit vector in index notation 
unit vector normal to a solid boundary 
above unit vector in index notation 
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N - (N x .N y .N z ) 


P 


q 

r 

R 


R min 


s 


S 



t 

•f\ 

t 


T 


u,v,w 


u i 


V 


x,y,z 


unit vector normal to the shock wave surface 
[expressed in the (x\y 1 ,z' )-sy$tem] 
vector normal to either a wave surface or a 
stream surface 
pressure 

velocity magnitude 
radial position of a point 
gas constant 
cowl lip radius 

distance from streamline base point to nearest 
point on convex hull 

either entropy per unit mass, or arc length 

temperature base in Sutherland's formula 

vector in the wave surface and normal to the 

bi characteristic direction 

time or time- like parameter 

unit vector along the space curve defined by 

the intersection of the shock wave with a 

meridional plane 

absolute temperature 

velocity components in the x, y, and z 

directions, respectively 

velocity in index notation 

velocity vector 

cartesian coordinates of base coordinate system 
base system coordinates in index notation 
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x\y\z' 


cartesian coordinates of local coordinate 


GREEK SYMBOLS 
a 


a V*i 


7 


5 


ij 


n 


0 


K 

A 

y 

t 

p 

o 

4> 

$ 


system 


either the angle of attack* or the angle sub- 

A 

tended by the unit vector Si and the z'-axis 
unit vectors used in the parameterization of 
the characteristic equations 
specific heat ratio 
Kronecker delta 

second coefficient of viscosity 
either the angle used in the parameterization 
of the characteristic equations* or the angle 
subtended by a meridian and the (x*y)-plane 
thermal conductivity 

term in the wave surface compatibility relation 
dynamic viscosity 
thermodynamic parameter 
density 

term in the noncharacteristic relation 

A 

angle subtended by the unit vector t and the 
x ' -axi s 

either the viscous dissipation function, or a 
term in the wave surface compatibility relation 
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SUBSCRIPTS 


waMiMmiwlBwiimri iir«tassafc **•*►«' 


\r ■■■■'* ' J ' 1 £ 1 


U.k 


x,y,z 


oo 


rectangular cartesian coordinate indices ranging 
from 1 to 3 

denotes either partial differentiation with 
respect to x, y, and z, or the x, y, and z 
components of a vector 
free-stream conditions 


OPERATORS 


D( )/Dt 

n 

n 


material derivative 

vector 

unit vector 
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